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METHOD AND SYSTEM FOR TRACKING 
MULTIPLE REGIONAL OBJECTS BY 
MULTI-DIMENSIONAL RELAXATION 

RELATED APPLICATIONS 

This application is a continuation-in-part of U.S. patent 
application Sen No. 08/404,024, filed Mar. 14. 1995 and 
issued Jul. 16, 1996 as U.S. Pat. No. 5,537,119, which is a 
continuation-in-part of U.S. patent application Ser. No. 
08/171 327 filed Dec. 21, 1993, now U.S. Pat. No. 5,406, 
289. 

FIELD OF THE INVENTION 

The invention relates generally to computerized tech- 
niques for processing data obtained from radiation reflec- 
tions used to track multiple discrete object. 

BACKGROUND OF THE INVENTION 

The invention relates generally to computerized tech- 
niques for processing data obtained from radar to track 
multiple discrete objects. 

There are many situations where the courses of multiple 
objects in a region must be tracked. Typically, radar is used 
to scan the region and generate discrete images or "snap- 
shots" based on sets of returns or observations. In some 
types of tracking systems, all the returns from any one object 
are represented in an image as a single point unrelated to the 
shape or size of the objects. 'Tracking" is the process of 
identifying a sequence of points from a respective sequence 
of the images that represents the motion of an object. The 
tracking problem is difficult when there are multiple closely 
spaced objects because the objects can change their speed 
and direction rapidly and move into and out of the line of 
sight for other objects. The problem is exacerbated because 
each set of returns may result from noise as well as echoes 
from the actual objects. The returns resulting from the noise 
are also called false positives. Likewise, the radar will not 
detect all echoes from the actual objects and this phenomena 
is called a false negative or "missed detect" error. For 
tracking airborne objects, a large distance between the radar 
and the objects diminishes the signal to noise ratio so the 
number of false positives and false negatives can be high. 
For robotic applications, the power of the radar is low and 
as a result, the signal to noise ratio can also be low and the 
number of false positives and false negatives high. 

In view of the proximity of the objects to one another, 
varied motion of the objects and false positives and false 
negatives, multiple sequential images should be analyzed 
collectively to obtain enough information to properly assign 
the points to the proper tracks. Naturally, the larger the 
number of images that are analyzed, the greater the amount 
of information that must be processed. 

While identifying the track of an object, a kinematic 
model describing the object's location, velocity and accel- 
eration may be generated. Such a model provides the means 
by which the object's future motion can be predicted. Based 
upon such a prediction, appropriate action may be initiated. 
For example, in a military application there is a need to track 
multiple enemy aircraft or missiles in a region to predict 
their objective, plan responses and intercept them. 
Alternatively, in a commercial air traffic control application 
there is a need to track multiple commercial aircraft around 
an airport to predict their future courses and avoid collision. 
Further, in these and other applications, such as robotic 
applications, may use radar, sonar, infrared or other object 
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detecting radiation bandwidths for tracking objects. In 
particular, in robotic applications reflected radiation can be 
used to track a single object which moves relative to the 
robot (or vice versa) so the robot can work on the object. 

5 Consider the very simple example of two objects being 
tracked and no false positives or false negatives. The radar, 
after scanning at time t,, reports objects at two locations in 
a first observation set. That is. it returns a set of two 
observations {o u , o J2 }. At time U it returns a similar set of 

10 two observations {o 2l , o 22 } from a second observation set. 
Suppose from prior processing that track data for two tracks 
T A and T 2 includes the locations at t 0 of two objects. Track 
Tj may be extended through the points in the two sets of 
observations in any -of four ways, as may track T 2 . The 

15 possible extensions of T, can be described as: {T,, o n , o 21 }. 
{T,, o u . o 22 }, \T t . o 12 , o 21 } and {Tj, o 12 , o 22 }. Tracks can 
likewise be extended from T 2 in four possible ways 
including, {T 2 , o 12 , o 21 }. FIG. 1 illustrates these five (out of 
eight) possible tracks (to simplify the problem for purposes 

20 of explanation) . The five track extensions are labeled h n , 
h 12 . h 13 , h 14 , and h 21 wherein h u is derived from {T^ o u , 
°2i K n i2 * s derived from {Tj, o u , o 22 }, h 13 is derived from 
{T x , o 12 , o 21 ), h l4 is derived from {Tj, o l2 , o 22 ), and h 21 is 
derived from {T 2 , o a , o 2i }. The problem of determining 

25 which such track extensions are the most likely or optimal 
is hereinafter known as the assignment problem. 

It is known from prior art to determine a figure of merit 
or cost for assigning each of the points in the images to a 
track. The figure of merit or cost is based on the likelihood 

30 that the point is actually part of the track. For example, the 
figure of merit or cost may be based on the distance from the 
point to an extrapolation of the track. FIG. 1 illustrates costs 
8 21 8 22 21 modeled target characteristics. The function to 
calculate the cost will normally incorporate detailed char- 

35 actenstics of the sensor, such as probability of measurement 
error, and track characteristics, such as likelihood of track 
maneuver. 

FIG. 2 illustrates a two by two by two matrix, c, that 

40 contains the costs for each point in relation to each possible 
track. The cost matrix is indexed along one axis by the track 
number, along another axis by the image number and along 
the third axis by a point number. Thus, each position in the 
cost matrix lists the cost for a unique combination of points 

45 and a track, one point from each image. FIG. 2 also 
illustrates a {0, 1} assignment matrix, z, which is defined 
with the same dimensions as the cost matrix. Setting a 
position in the assignment matrix to "one" means that the 
equivalent position in the cost matrix is selected into the 

50 solution. The illustrated solution matrix selects the {h I4 . 
h 2 J solution previously described. Note that for the above 
example of two tracks and two snapshots, the resulting cost 
and assignment matrices are three dimensional. As used in 
this patent application, the term "dimension" means the 

5S number of axes in the cost or assignment matrix while size 
refers to the number of elements along a typical axis. The 
costs and assignments have been grouped in matrices to 
facilitate computation. 
A solution to the assignment problem satisfies two 

60 constraints — first, the sum of the associated costs for assign- 
ing points to a track extension is minimized and, second, if 
no false positives or false negatives exist, then each point is 
assigned to one and only one track. 
When false positives exist, however, additional hypotheti- 

65 cal track extensions incorporating the false positives will be 
generated. Further note that the random locations of false 
positives will, in general, not fit well with true data and such 
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additional hypothetical track extensions will result in higher 
costs. Also note that when false negative errors exist then 
the size of the cost matrix must grow to include hypothetical 
track extensions formulated with "gaps" (i.e.. data omissions 
where there should be legitimate observation data) for the 
false negatives. Thus, the second criteria must be weakened 
to reflect false positives not being as signed and also to 
permit the gap filler to be multiply assigned. With hypo- 
thetical cost calculated in this manner then the foregoing 
criteria for minimization will tend to materialize the false 
negatives and avoid the false positives. 

For a 3-dimensional problem, as is illustrated in FIG. 1, 
but with N x (initial) tracks. N 2 observations in scan 1. N 3 
observations in scan 2. false positives and negatives 
assumed, the assignment problem can be formulated as: 

"l A* *3 11.0] 

(z) Minimize: £ £ £ ^y^i'so 

(b) Subject to: £ £ W3 = *' = 1 "" Nl 

<2*1<3 B 1 
'I-l '3-1 

(d) 2]Z*lW S1, 1*3 = 1, - ^3 

<e) U&ty €{0,1} V;,^ 

where "c** is the cost and "z" is a point or observation 
assignment, as in FIG. 2. 

The minimization equation or equivalently objective 
function [1.0] (a) specifies the sum of the element by 
element product of the c and z matrices. The summation 
includes hypothesis representations w * m observation 
number zero being the gap filler observation. Equation 
[1.01] (b) requires that each of the tracks T x . . . . .T N be 
extended by one and only one hypothesis. Equation [1.0] (c) 
relates to each point or observation in the first observation 
set and requires that each such observation, except the gap 
filler, can only associate with one track but because of the 
"less than" condition it might not associate with any track. 
Equation [1.01] (d) is like [1.0] (c) except that it is appli- 
cable to the second observation set. Equation [1.0] (e) 
requires that elements of the solution matrix z be limited to 
the zero and one values. 

The only known method to solve Problem Formulation 
[1.0] exactly is a method called "Branch and Bound." This 
method provides a systematic ordering of the potential 
solutions so that solutions with a same partial solution are 
accessible via a branch of a tree describing all possible 
solutions whereby the cost of unexamined solutions on a 
branch are incrementally developed as the cost for other 
solutions on the branch are determined. When the develop- 
ing cost grows to exceed the previously known minimal cost 
(Le.. the bound) then enumeration of the tree branch termi- 
nates. Evaluation continues with a new branch. If evaluation 
of the cost of a particular branch completes, then that branch 
has lower cost than the previous bound so the new cost 
replaces the old bound. When all possible branches are 
evaluated or eliminated then the branch mat had resulted in 
the last used bound is the solution. If we assume that 
NjsN^N^n and that branches typically evaluate to half 
there full length, then workload associated with "branch and 
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bound** is proportional to (nllV:!) 2 . This workload is 
unsuited to real time evaluation. 

The Branch and Bound algorithm has been used in past 
research on the Traveling Salesman Problem. Messrs. Held 

5 and Karp showed that if the set of constraints was relaxed by 
a method of Lagrangian Multipliers (described in more 
detail below) then tight lower bounds could be developed in 
advance of enumerating any branch of the potential solution. 
By initiating the branch and bound algorithm with such a 

10 tight lower bound, significant performance improvements 
result in that branches will typically evaluate to less than half 
their full length. 

Messrs. Frieze and Yadagar in dealing with a problem 
related to scheduling combinations of resources, as in job. 

15 worker and work site, showed that Problem Formulation 
[1.0] applied. They further described a solution method 
based upon an extension of the Lagrangian Relaxation 
previously mentioned. The two critical extensions provided 
by Messrs. Frieze and Yadagar were: (1) an iterative proce- 

20 dure that permitted the lower bound on the solution to be 
improved (by "hill climbing" described below) and (2) the 
recognition that when the lower bound of the relaxed 
problem was maximized, then there existed a method to 
recover the solution of the non-relaxed problem in most 

25 cases using parameters determined at the maximum. The 
procedures attributed to Messrs. Frieze and Yadagar are only 
applicable to the 3-dimensional problem posed by Problem 
Formulation f 1.0] and where the cost matrix is fully popu- 
lated. However, tracking multiple airborne objects usually 

30 requires solution of a much higher dimensional problem 
FIGS. 1 and 2 illustrate an example where "look ahead** 
data from the second image improved the assignment accu- 
racy for the first image. Without the look ahead, and based 
only upon a simple nearest neighbor approach, the as sign - 

35 ments in the first set would have been reversed. Problem 
Formulation (1.0] and the prior art only permit looking 
ahead one image. In the prior art it was known that the 
accuracy of assignments will improve if the process looks 
further ahead, however no practical method to optimally 

40 incorporate look ahead data existed. Many real radar track- 
ing problems involve hundreds of tracks, thousands of 
observations per observation set and matrices with dimen- 
sions in the range of 3 to 25 including many images of look 
ahead. 

45 It was also known that the data assignment problem may 
be simplified (without reducing the dimension of the assign- 
ment problem) by eliminating from consideration for each 
track those points which, after considering estimated limits 
of speed and turning ability of the objects, could not 

50 physically be part of the track. One such technique, denoted 
hereinafter the "cone method." defines a cone as a continu- 
ation of each previously determined track with the apex of 
the cone at the end of the previously defined track. The 
length of the cone is based on the estimated maximum speed 

55 of the object and the size of the arc of the cone is based on 
the estimated maximum turning ability of the object Thus, 
the cone defines a region outside of which no point could 
physically be part of the respective track. For any such 
points outside of the cones, an infinite number could be put 

60 in the cost matrix and a zero could be preassigned in the 
assignment matrix. It was known for the tracking problem 
that these elements will be very common in the cost and 
selection matrices (so these matrices are "sparse**). 
Il was also known in the prior art that one or more tracks 

65 which are substantially separated geographically from the 
other tracks can be separated also in the assignment prob- 
lem. This is done by examining the distances from each 
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point to the various possible tracks. If the distances from one observations in a particular scan to which the constraint 

set of points are reasonably short only in relation to one set is related. In general, there is a constraint for each 

track, then they are assigned to that track and not further observation of the scan, wherein the constraint indi- 

considered with the remainder of the points. Similarly, if a cates the number of track extensions to which the 

larger group of points can only be assigned to a few tracks, 5 observation may belong. 

then the group is considered a different assignment problem. In particular, the first optimization problem is formulated. 

Because the complexity of assignment problems increases generated or defined as an M -dimensional assignment prob- 

dramatically with the number of possible tracks and the total iem wherein there are M constraint sets in the first collection 

number of points in each matrix, this partitioning of the of constraint sets (i.e.. there are M scans being examined) 

group of points into a separate assignment problem and 10 and the first objective function minimizes a total cost for 

removal of these points from the matrices for the remaining assigning observations to various track extensions wherein 

points, substantially reduces the complexity of the overall terms are included in the cost, such that the terms have the 

assignment problem. figures of merit or costs for hypothesized combinations of 

A previously known Multiple Hypothesis Testing (MHT) assignments of the observations to the tracks. Subsequently, 
algorithm (see Blackman, Multiple-Target Tracking with is the formulated M-dimensional assignment problem is 
Radar Applications, Chapter 10, Artech House Norwood solved by reducing the complexity of the problem by 
MA, 1986) related to formulation and scoring. The MHT generating one or more optimization problems each having 
procedure describes how to formulate the sparse set of all a lower dimension and then solving each lower dimension 
reasonable extension hypothesis (for FIG. 1 the set ^h u . . . optimization problem. That is, the M-dimensional assign- 
h 24 }) and how to calculate a cost of the hypothesis {T,. o ly , 20 ment problem is solved by solving a plurality of optimiza- 
o^j based upon the previously calculated cost for hypoth- tion problems each having a lower number of constraint sets, 
esis {T,. o^}. The experience with the MHT algorithm. The reduction of the M-dimensional assignment problem 
known in the prior art, is the basis for the assertion that look to a lower dimensioned problems is accomplished by relax- 
ahead through k sets of observations results in improved ing the constraints on the points of one or more scans 
assignment of observations from the first set to the track. 25 thereby permitting these points to be assigned to more than 

In theory, the MHT procedure uses the extendable costing one track extension. In relaxing the constraints, terms having 

procedure to defer assignment decision until the accumu- penalty factors are added into the objective function thereby 

lated evidence supporting the assignment becomes over- increasing the total cost of an assignment when one or more 

whelming. When it makes the assignment decision it then points are assigned to more than one track. Thus, the 

eliminates all potential assignments invalidated by the deci- 30 reduction in complexity by this relaxation process is itera- 

sion in a process called **pruning the tree." In practice, the tively repeated until a sufficiently low dimension is attained 

MHT hypothesis tree is limited to a fixed number of gen- such that the lower dimensional problem may be solved 

erations and the overwhelming evidence rule is replaced by directly by known techniques. 

a most likely and feasible rule. This rule considers each track In one embodiment of the invention, each k-dimensional 

independently of others and is therefore a local decision rule. 35 assignment problem 2<k2lM, is iteratively reduced to a k-1 

A general object of the present invention is to provide an dimensional problem until a 2 -dimensional problem is speci- 

efficient and accurate process for assigning each point object fied or formulated. Subsequently, the 2 -dimensional prob- 

in a region from multiple images to a proper track and then lem formulated is solved directly and a "recovery" technique 

taking an action based upon the assignments. is used to iteratively recover an optimal or near-optimal 

A more specific object of the present invention is to 40 solution to each k-dimensional problem from a derived 

provide a technique of the foregoing type which determines (k-1) dimensional problem k=2, 3, 4, . . .M. 

the solution of a k-dimensional assignment problem where In performing each recovery step (of obtaining a solution 

"k" is greater than or equal to three. to a k-dimensional problem using a solution to a (k-1)- 

SUMMARY OF THE INVENTION dimensional problem) an auxiliary function, is utilized. In 

45 particular, to recover an optimal or near-optimal solution to 

The present invention relates to a method and apparatus a ^.dimensional problem, an auxiliary function. HV*. 

for tracking objects. In particular, the present invention 5 M is specified and a region or domain is determined 

tracks movement or trajectories of objects by analyzing wherein this function is maximized, whereby values of the 

radiation reflected from the objects, the invention being region determine the penalty factors of the (k-1)- 

especially useful for real-time tracking in noisy environ- 50 dimensional problem such that another 2-dimensional prob- 

ments. lem can be formulated which determines a solution to the 

In providing such a tracking capability, a region contain- k-dimensional problem using the penalty factors of the 

ing the objects is repeatedly scanned to generate a multi- (k-l)-dimensional problem. 

plicity of sequential images or data observation sets of the Each Auxiliary function 4\ of both lower dimensional 

region. One or more points (or equivalently observations), in 55 problem penalty factors and a solution at the dimension k at 

each of the images or observation sets are detected wherein which the penalized cost function is solved directly 

each such observation either corresponds to an actual loca- (typically a 2-dimensional problem). Further, in 

tion of an object or is an erroneous data point due to noise. determining, for auxiliary function determined, and 

Subsequently, for each observation detected, figures of merit utilized to identify the peak region. Thus, gradients are used 

or costs are determined for assigning the observation to each $o for each of the approximation functions 4^ 4* 4 determining 

of a plurality of previously determined tracks. * Afterwards, penalty factors penalty factors until tm-1 is used in deter- 

a first optimization problem is specified which includes: mining the penalty factors for the M-l dimensional prob- 

(a) a first objective function for relating the above men- lem. Subsequently, once the M-dimensional problem is 
tioned costs to potential, track extensions through the solved (using a 2-dimensional problem to go from an (M-l) 
detected observations (or simply observations); and 65 dimensional solution to an M-dimensional solution), one or 

(b) a first collection of constraint sets wherein each more of the following actions are taken based on the track 
constraint set includes constraints to be satisfied by the assignments: sending a warning to aircraft or a ground or sea 
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facility, controlling air traffic, controlling anti-aircraft or Other features and benefits of the present invention will 

anti-missile equipment, taking evasive action, working on become apparent from the detailed description with the 

one of the objects. accompanying drawings contained hereinafter. 

According to one feature of this first embodiment of the 

present invention, the following steps are also performed 5 BRIEF DESCRIPTION OF THE FIGURES 

before the step of defining the auxiliary function. A prelimi- nG j i$ a of images or ^ sets generated by a 

nary auxiliary function is defined for each of the lower scan of a region afld possible n&tks within me images or 

dimensional problems having a dimension equal or one ^ seU according t0 tne pr ior art 

greater than the dimension at which the penalized cost ^ _ . # _ „ ^ 

r . . , * ■ • t ™. i • •!■ a FIG. 2 illustrates cost and assignment matrices for the 

function is solved directly. The preliminary auxiliary func- 10 * * T~!" - trt * . r Qrt 

^- r i \, M „ 0 ii, r 0 „h o data sets of FIG. 1 according to the prior art 

tion is a function of lower order penalty values and a fr r 

solution at the dimension at which the penalized cost func- FIG. 3 is a block diagram of the present invention, 
tion was solved directly. In deterrnining a gradient of the FIG. 4 is a flow chart of a process according to the prior 
preliminary auxiliary function, step in the direction of the art for solving a 3-dimensional assignment problem, 
gradient to identify a peak region of the preliminary auxil- 15 FIG. 5 Is a flow chart of a process according to the present 
iary function and determine penalty factors at the peak invention for solving a k-dimensional assignment problem 
region. Iterativeiy repeat the defining, gradient determining. where "k" is greater than or equal to 3. 
stepping and peak determining steps to define auxiliary nQ 6 h a h of various functions used t0 explain foe 
functions at successively higher dimensions until the auxil- pre sent invention, 
iary function at 6-18 (k-1) dimension is determined In an 2 o ™^ *, • / ^ , j- r 
Snative second emblem of the present invention. f ™ ; 7 15 . Wock : dm f m ^tuZ^Z 
instead of reducing the dimentiality of the M-dimensional *" ^ the k-dimensional assignment problem where 
assignment problem by a single dimension at a time, a k 15 or e 1 ual t0 3 ' „ , , . 
plurality of dimensions are relaxed simultaneously. This new FIG. 8 is a flowchart describing the procedure for solving 
strategy has the advantage that when the M-dimensional 25 a n-dimensional assignment problem according to the sec- 
problem is relaxed directly to a 2-dimensional assignment ond embodiment of the invention, 
problem, then all computations may be performed precisely DETAILED DESCRIPTION OF THE 
without utilizing an auxiliary function such as ¥ k embodi- PREFERRED EMBODIMENTS 
ment More particularly, the second embodiment solves the 

first optimization problem (i.e.. the M-dimensional assign- 30 Referring now to the other figures in detail wherein like 

ment problem) by specifying (i.e.. creating, generating. reference numerals indicate like elements throughout the 

formulating and/or defining) a second optimization problem. several views. FIG. 3 illustrates a system generally desig- 

The second optimization problem includes a second objec- nated 100 for implementing the present invention. System 

rive function and a second collection of constraint sets 100 comprises, for example, a radar station 102 (note sonar, 

wherein: 35 microwave, infrared and other radiation bandwidths are also 

a) the second objective function is a combination of the contemplated) for scanning a region which may be. for 
first objective function and penalty factors or terms example, an aerial region (in aerial surveillance 
determined for incorporating M-m constraint sets of the applications) or a work region (in robotic applications) and 
first optimization problem into the second objective generating signals indicating locations of objects within the 
function* 40 re gt° n - The signals are input to a converter 104 which 

b) the constraint sets of the second collection include only converts * e si S nals to P° ints or observations in which 
m of the constraint sets of the first collection of each object (or false positive) is represented by a single 
constraints P° mt 0Ut P ut of toe conveiter & in P ut t0 readable by 

wherein 2£m£M-l. Note that, once the second optimiza- a computer 106. As described in more detail below, the 
tion problem has been specified or formulated, an optimal or 45 computer 106 assigns the points to respective tracks and 
near-optimal solution is determined and that solution is used then displays the tracks and extrapolations ; of the tracks on 
in specifying (i.e.. creating, generating, formulating and/or a monitor 110. Also, the computer 106 determines an 
defining) a third optimization problem of M-m dimensions appropriate action to take based on the tracks and track 
(or equivalently constraint sets). The third optimization extensions. For example, in a commeraal application at an 
problem is subsequently solved by decomposing it using the 50 airport, the computer can determine if two aircraft being 
same procedure of this second embodiment as was used to tracked are on a collision course and if so signal a trans- 
decompose the first optimization problem above. Thus, a nutter 112 to warn each aircraft, or if a scheduled take-off 
plurality of instantiations of the third optirnization problem w ^ P° se the risk of collision, delay the take-off. For a 
are specified, each successive instantiation having a lower rnilitary application on a ship or base, the computer can 
number of dimensions, until an instance of the third opti- 55 determine subsequent coordinates of enemy aircraft and 
mization problem is a two dimensional assignment problem send the coordinates to an anUaircraft gun or missile 120 via 
which can be solved directly. Subsequently, whenever an a communication channel 122. In a robotic apphcauon. the 
instance of the third optimization problem is solved, the computer controls the robot to work on the proper object or 
solution is used to recover a solution to the instance of the portion of the object 

first optimization problem from which this instance of the 60 The invention generates k-dimensional matrices where k 

third optimization was derived. Thus, an optimal or near- is the number of images or sets of observation data in the 

optimal solution to the original first optimization problem look ahead window plus one. Then, the invention formulates 

may be recovered through iteration of the above steps. a k-dimensional assignment problem as in [1.0] above. 

As mentioned above, the second embodiment of the The k-dimensional assignment problem is subsequently 

present invention is especially advantageous when m=2, 65 relaxed to a (k-l)-dimensional problem by incorporating 

since in this case all computations may be performed one set of constraints into the objective function using a 

precisely and without utilizing auxiliary functions. Lagrangian relaxation of this set Given a solution of the 
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(k-l)-dimensionaI problem, a feasible solution of the observation with k-1 gap fillers, e.g., Zq ... 0 <,°-- '°.i p *0- 
k-dimensional problem is then reconstructed. The (k-1)- Thus, the resulting generalization of Problem Formu- 
dimensional problem is solved in a similar manner, and the lation 1 1.0| without the "less than" complication within 
process is repeated until it reaches the two-dimensional the constraints is the following k-Dimensional Assign- 
problem that can be solved exactly. The ideas behind the 5 ment Problems in which k£3: 
Lagrangian relaxation scheme are outlined next. 
Consider the integer programming problem . % & 11,11 

. /x r Minimize: 2, - L e, i~'i°i 

Minimize v(z) =c z i t -o 

Subject to: Az ^b io 
Bz ^d 

z, is an integer for i el 
11.0.11 

where the partitioning of the constraints is natural in some * t "j-i " >+1 N k 
sense. Given a multiplier vector ulO, the Lagrangian relax- £ ... ]T X "2 El- 
ation of 1 1.0.1] relative to the constraints Bz^d is defined to '>-i°° 'j+i* 0 
be 

(DOTT 

Subject to: Az^b 

z ( is an integer for i el 20 Wj Nk l 

|1:0.2| £ ^ = i. 

If the constraint set Bz^d is replaced by Bz=d. the nonne- 'i=*> 

gativiry constraint on u is removed. M =C r z+u r (Bz-d) is the * t <= (0, l|, Vi m n=\ k 

Lagrangian relative to the constraints Bz=d, and hence the ^ 

name Lagrangian relaxation. The following fact gives the . 

relationship between the objective functions of the original where c and z ™ similarly dimensioned matrices represent- 

and relaxed problems in S costs and hypothetical assignments. Note, in general, for 

FACT A.l. If z is'an optimal soluUon to | 1.0.1], then tracing problems these matrices are sparse. 

. ^ ,\ . „ ^ A Tf y . . ' Afl1Aor After formuIaUng the k-dimensional assignment problem 

^v(z) for all I uJC If an opUnial solution z of 11.02] is ^ solyes th | resultin P ^ 

feasible for f 1.0.2], then z is an optimal solution for [1.0.1] ^ ^ J Xq ^ outputs required by 

ana <PAlgor lthm 112 U2 and m For ttdl observation set 0| received from 

converter 104 at time t ( where i=l , . . , Nj, the computer 

k k = o processes O, in a batch together with the other observation 

35 sets 0^ +l O, and the track T,.*, i.e.. Tj is the set of all 

_ c . . - tracks that have been defined up to but not including (X. 

converging to the solution u of Maximize { 0^0} and a (Note Md designations refer t0 me vector of dements 

corresponding sequence of feasible solutions {Z.^Oof for me mdicated time i e ^ me set of all observations in the 

(1.0.1] as follows: scan 0f Uacks cx i st i„g at tne tf mt% e tc.) The result of this 

1. Generate in initial approximation 40 processing is the new set of tracks T^ +l and a set of cost 

2. Given u A , choose a search direction s k and a search weighted possible solutions indicating how the tracks might 
distance a k H* k a* k *Y k extend to the current time t ( . At time t i+J the batch process 

estimate u k by u^^u^-hx* k is repeated using the newest observation set and deleting the 

3. Given u A+1 and a feasible solution t k+l (u A+A ) of [1.0.2], oldest. Thus, there is a moving window of observation sets 
recover a feasible solution z t+1 (u t+1 ) of the integer 45 which is shifted forward to always incorporate the most 
programming problem 1 1 0.1 1 recent observation set. The effect is that input observation 

4. Check the termination criteria. If the algorithm is not sets ^ e ' eused for c y c,es and * en ° K n observation 
finished, set k=k+l and return to Step 2. Otherwise, set s k " th re " se each observation within the observation set 
» . , A ..^ is integrated into a track. 

terminate the algorithm. , A _ ... _ A . . . 

Tf - . I .*. aPM An M 50 FIG. 7 illustrates various processes implemented upon 

If z is an optimal solution of 1.0.1], then . 4 _ . _ . . * .... *T 

_ receipt of each observation set. Except for the addition of the 

^*-^- v ( z )-( z a) , r*^«i- .1 k-dimensional assignment solving process 300 and the 

Since the optimal soluuon z=2(u) of [ 1.0.2] is usuaUy not a modification t0 scoring process 154 to buiW ^ structures 

feasible solution of [1.0.1] for any choice of the multipliers, suitaWe fof pgQcm 300 ^ processcs m nG 7 m 5ased 

ss upon prior art. The following processes 150 and 152 extend 

* previously defined tracks h t _ lm based on new observations. 

* a a Gate formulation and output process 156 determines, for 

each of the previously defined tracks, a zone wherein the 
track may potentially extend based on limits of velocity, 

« 60 maneuverability and radar precision. One such technique to 

a a accomplish this is the cone method described previously. 

The definition of the zone is passed to gating process 150. 
When a new observation set O i is received, the gating 

p^p^k in z*, i* in problem [1.1] below indicates that process 150 will match each member observation with the 

the track extension represented by z*^ i k includes a 65 zone for each member of the hypothetical set h M . After all 

false positive in the p"' observation set Note that this input observations from 0 ( are processed the new hypothesis 

implies that a hypothesis be formed incorporating an set h ( is generated by extending each track of the prior set of 
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hypothetical tracks h, x either with missed detect gap fillers observation identifiers, and passes them to k-dimensional 
or with'all new observation elements satisfying the track's assignment problem solving process 300. 
zone. This is a many to many matching in that each The assignment solving process 300 is described below, 
hypothesis member can be extended to many new observa- Its output is simply the list of assignments which constitute 
tions and each new observation can be used to extend many 5 the most likely solution of the problem described by Equa- 
hypotheses. It however, is not a full matching in that any tion [1.1]. Note that both gate formulation and output 
hypothesis will neither be matched to all observations nor processes 156 and 162 use (at different times) the list of 
vice versa. It is this matching characteristic that leads to the assignments to generate the updated track history T, to 
sparse matrices involved in the tracking process. eliminate or prune alternative previous hypotheses that are 
Subsequently, gating 150 forwards the new hypothesis set h f 10 prohibited by the actual assignments in the list, and subse- 
to filtering process 152. Filtering process 152 determines a ' quently to output any required data. Also note that when one 
smooth curve for each member of h ( . Such a smooth curve of the gate formulation and output processes 156 and 162 
is more likely than a sharp turn from each point straight to accomplish these tasks, the process will subsequently gen- 
the next point. Further, the filtering process 152 removes erate and forward the new set of gates for each remaining 
small errors that may occur in generating observations. Note 15 hypothesis and the processes will then be prepared to receive 
that in performing these tasks, the filtering process 152 the next set of observations. In one embodiment, the loop 
preferably utilizes a minimizauon of a least squares test of described here will generate zones for a delayed set of 
the points in a track hypothesis or a Kalman Filtering observations rather than the subsequent set This permits 
approach. processes 156 and 162 to operate on even observation sets 
As noted above, the foregoing track extension process 20 while the scoring step 154 and k-dimensional solving pro- 
requires knowledge of a previous track. For the initial cess operate on odd sets of observations, or vice versa, 
observations, the following gating process 158 and filtering The assignment solving process 300 permits the present 
process 160 determine the "previous track" based on initial invention to operate with window size of k-1 for k!3. The 
observations. In deterrnining the initial tracks, the points upper limit on k depends only upon the computational power 
from the first observation set form the beginning points of all 25 of the computer 106 and the response time constraints of 
possible tracks. After observation data from the next obser- system 100. The k-1 observation sets within the processing 
vation set is received, sets of simple two point straight line window plus the prior track history result in a k-dimensional 
tracks are defined. Then, promotion, gate formulation, and Assignment Problem as described by Problem Formulation 
output step 162 determines a zone in which future exten- [1.1]. The present invention solves this generalized problem 
sions are possible. Note that filtering step 160 uses curve 30 including the processes required to consider false positives 
fitting techniques to smooth the track extensions depending and negatives, and also the processes required to consider 
upon the number of prior observations that have accumu- sparse matrix problem formulations, 
lated in each hypothesis. Further note that promotion, gate L A FIRST EMBODIMENT OF THE k-DIMENSIONAL 
formulation and output process 162 also determines when ASSIGNMENT SOLVER 300 

sufficient observations have accumulated to form a basis for 35 In describing a first embodiment of the k-dimensional 
promoting the track to processes 150 and 152 as described assignment solver 300. it is worthwhile to also discuss the 
above process of FIG. 4 which is used by the solver 300. FIG. 4 
The output of each of the filtering processes 152 and 160 illustrates use of the Frieze and Yadagar process as shown in 
is a set of hypothetical track extensions. Each such extension prior art for transforming a 3-dimensional assignment prob- 
contains the hypothetical initial conditions (from the previ- 40 lem into a 2-dimensional assignment problem and then use 
ous track), the list of observations incorporated in the a hill climbing algorithm to solve the 3-dimensional assign- 
extension, and distance between each observation and the ment problem. The solver 300 uses a Lagrangian Relaxation 
filtered track curve. Scoring process 154 determines the technique (well known in the art) to reduce the dimension of 
figure of merit or cost of an observation being an extension an original k-dimensional assignment problem (k>3) down 
of a track. In one embodiment the cost is based on the 45 to a 3-dimensional problem and then use the process of FIG. 
above-mentioned, distance although the particular formula 4 to solve the 3-dimensional problem Further note that the 
for determining the cost is not critical to the present inven- Lagrangian Relaxation technique is also utilized by the 
tion. A preferred formula for determining the cost utilizes a process of FIG. 4 and that in using this technique the 
negative log likelihood function in which the cost is the requirement that each point is assigned to one and only one 
negative of the sum of: (a) the logs of the distances nor- 50 track is relaxed. Instead, an additional cost, which is equal 
malized by sensor standard deviation parameters, and (b) the to a respective Lagrangian Coefficient "u". is added to the 
log likelihoods for events related to: track initiation, track cost or objective function 1 1.0] (a) whenever a point is 
termination, track maneuver, false negatives and false posi- assigned to more than one track. This additional cost can be 
tives. Note that track maneuvers are detected by comparing picked to weight the significance of each constraint violation 
the previous track curve with the current extension. Further 55 differently, so this additional cost is represented as a vector 
note that some of the other events related to. for example, of coefficients u which are correlated with respective obser- 
false negatives and false positives are detected by analyzing vation points. Hill climbing will then develop a sequence of 

the relative relationship of gap fillers in the hypothesis. Lagrangian Coefficients sets designated (u 0 ">* u /*i 

Thus, after determining that one of these events occurred, a u p ). That correspond to an optimum solution of the 

cost for it can be determined based upon suitable statistics 60 2-dimensional assignment problem The assignments at this 

tables and system input parameters. The negative log like- optimum solution are then used to "recover" the assignment 

lihood function is desirable because it permits effective solution of the 3-dimensional assignment problem 

integration of the useful components. Copies of the set of In step 200 of FIG. 4. initial values are selected for the u 0 

hypothetical track extensions which are scored are subse- coefficients. Because the Lagrangian Relaxation process is 

quently passed directly to one of the gate formulation and 65 iterative, the initial values are not critical and are all initially 

output steps 156 and 162. Note that the scoring process 154 selected as zero. In step 202, these additional costs are 

also arranges the actual scores in a sparse matrix based upon applied to the objective function [ 1.0] (a). With the addition 
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of the costs 4 V\ the goal is still to assign the points which 
minimize the total cost. This transforms Equation 1 1.01 (a), 
written for k=3 and altered to exclude mechanisms related to 
false positives and negatives, into objective function |2.1| 
(a) In the first iteration it is not necessary to consider the "u" 
matrix because all values are set to zero. To relax the 
requirement that each point be assigned to one and only one 
track, the constraint Equation 1 1.0] (d) is deleted, thereby 
permitting points from the last image to be assigned to more 
than one track- Note that while any axis can be chosen for 
relaxation, observation constraints are preferably relaxed. 
The effect of this relaxation is to create a new problem which 
must have the same solution in the first two axes but which 
can have a differing solution in the third axis. The result is 
constraints |2.1| (b-d). 



(a) Minimize: 



Z Z K c 'i'2<3 " % H*l*3 



12.11 



(b) Subjccl to: £ £ z*,^ ■■ 



(d) 



ZZ*^* 1 - 

*,V3 €{0, U 



i 



k = 1 N 2 



Step 204 then generates from the 3 -dimensional problem 
described by Problem Formulation 11.0] a new 
2-dimensional problem formulation which will have the 
same solution for the first two indices. As Problem Formu- 
lation 1 2.1 1 has no constraints on the 3^ axis, any value 
within a particular 3^ axis can be used in a solution, but 
using anything other than the minimum value from any 3-rd 
axis has the effect of increasing solution cost. Conceptually, 
the effect of step 204 is to change the 3-dimensional arrays 
in Problem Formulation |2.11 into 2-dimensional arrays as 
shown in Problem Formulation |2.2| and to generate the new 
2-dimensional matrix m ( J 2 defined as shown in Equation 
12.3]. 



(a) Minimize: 



tt 



Min: 



12.21 



Vij 



(b) Subject to: 



zz^* 



1, I'i = 1. 



(c) z Z*^ sl * h ~ l Ni 

(d) Zi^ € |0, 1} V Z4 il2 

m ili2 « Min: arg minimize (c^ - u i} \ i = 1, . 
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<l> u 



10 



15 



20 



25 



30 



coefficients u^j whose corresponding Problem Formulation 
1 2.2 1 problem solution is a cost value closer to the peak of 

dimensional assignment problem requires solving the Equa- 
tion 2.2 problem corresponding to u p . 

In step 206. the two dimensional problem is solved 
directly using a technique known to those skilled in the art 
such as Reverse Auction for the corresponding cost and 
solution values. This is the evaluation of one point on the 
surface or for the first iteration 4>U 0 

Thus, after this first iteration, the points have been 
assigned based on all "u" values being arbitrarily set to zero. 
Because the "u" values have been arbitrarily assigned, it is 
unlikely that these assignments are correct and it is likely 
that further iterations are required to properly assign the 
points. Step 208 determines whether the points have been 
properly assigned after the first iteration by determining if 
for this set of assignments whether a different set of "u" 
values could result in a higher total cost. Thus, step 208 is 
implemented by determining the gradient of objective func- 
tion |2.2 1 (a) with respect to u ; . If the gradient is substan- 
tially non-zero (greater than a predetermined limit) then the 
assignments are not at or near enough to the peak of the d>u 



35 



7+i 



45 



12.31 



The cost or objective function for the reduced problem as 
defined by [2.21 (a), if evaluated at all possible values of u 
is a surface over the domain of U 3 . This surface is referred 
to as <tHi 



Hill climbing Step 212 determines the u J+l values based 
upon the u y values, the direction resulting from protecting 
the previous gradient into the U 3 domain, and a step size. 
40 The solution value of the 2-dimensional problem is the set 
of coefficients that minimize the 2-dimensional problem and 
the actual cost at the minimum. Those coefficients aug- 
mented by the coefficients stored in m,^ permit the evalu- 
ation (but not the minimization) of the cost term in Problem 
Formulation 1 2. 1 1. These two cost terms are lower and upper 
bounds on the actual minimized cost of the 3-dimensional 
problem, and the difference between them in combination 
with the gradient is used to compute the step size. 

With this new set of "u" values, steps 202-210 are 
repeated as a second iteration. Steps 212 and 202-210 are 
repeated until the gradient as a function of u determined in 
step 208 is less than the predetermined limit. This indicates 
that the up values which locate the peak area of the <1> u 
surface are determined and that the corresponding Problem 
Formulation [2.2] has been solved. Step 214 will attempt to 
use the assignments that resulted from this particular 
2-dimensional assignment problem to recover the solution of 
the 3-dimensional assignment problem as described below. 
If the limit was chosen properly so that the "u M values are 
close enough to the peak, this recovery will yield the proper 
set of assignments that rigidly satisfies the constraint that 
each point be assigned to one and only one track. However, 
if the 44 u" values are not close enough to the peak, then the 
limit value for decision 210 is reduced and the repetition of 
steps 212 and 202-210 is continued. 



50 



55 



60 



65 





and / = 2. ...A - 1 



(d> Z-I^'C,- 1 - 
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Step 214 recovers the 3-dimensional assignment solution 
by using the assignment values determined on the last ^ ^ I3.1J 

iteration through step 208. Consider the 2-dimensional z < a > Minimize: v t tf)= l/-^...^,...^ 
assignment matrix to have Ts in the locations specified by 

the list L,=(a,. b^„. If the 3-dimensional z matrix is 5 (b) ^ % A > . = t ^ 

specified by placing l's at the location indicated by the list £j £4 * 

LjKa,. b f . ra^)^ then the result is a solution of Problem ^ ^ ^ 

Formulation [2.1]. Let L 3 =(m fl * I .) Afel be the list formed by (c) Z " Z Z "Z^i 'i = 1 

the third index. If each member of L 3 is unique then the L 2 10 m- 0 v-j- 0 

solution satisfies the third constraint so it is a solution to for /■ = 1, ... at> 

Problem Formulation 1 1.0]. When this is not the case, 
recovery determines the minimal substitutions required 
within list L 3 so that it plus L : will be a feasible solution. i.e., 

a solution which satisfies the constraints of a problem 15 ( — 0 
formulation, but which may not optimize the objective . m . . 

(c) 2J. € |0, 1J, for all fi... li 

function of the problem formulation. This stage of the l "" * 

recovery process is formulated as a 2-dimensional assign- 
ment problem: Form a new cost matrix [c^]"^ where 20 and where ^ is ^ cost matrix [c*. ... .J which is a function 
tifVofaj for j=l • • • N / me N * term is me 10X31 number of ^ di^ce between the observed point z* and the 
of cost elements in the selected row of the 3-dimensional smoothed track determined by the filtering step, and v k is the 
cost matrix. Attempt to solve this 2-dimensional problem for cost f unct j 0IL This set of equations is similar to the set 
the minimum using two constraints sets. If a feasible solu- presented in Problem Formulation [1.1] except that it 
tion is found then the result will have the same form as list 25 indudes me su b SC ript and superscript k notation. Thus, in 
L x . Replace the first set of indexes by the indicated (a ( . b ; ) solving the M-dimensional Assignment Problem the inven- 
pairs taken from list L : and the result will be a feasible tion reduces this pro blem to an M-l dimensional Assign- 
solution of Problem Formulation [1.1]. If no feasible solu- ment p^oiem and then to an N-2 dimensional Assignment 

tion to the new 2-dimensional problem exists then further 30 problem, etc. Further, the symbol ke{3 M} customizes 

effort to locate the peak of O u Problem Formulation [3.1] to a particular relaxation level. 

That is, the notation is used to reference data from levels 

* relatively removed as in c** 1 are the cost coefficients which 

existed prior to this level of relaxed coefficients c*. Note that 
actual observations are numbered from 1 to N, where N, is 
the number of observations in observation set i. Further note 
that the added 0 observation in each set of observations is the 
unconstrained "gap filler." This element serves as a filler in 
* 40 substituting for missed detects, and a sequence of these 

* elements including only one true observation represents the 

possibility that the observation is a false positive. Also note 
that by being unconstrained a gap filler may be used in as 
many hypotheses as required. 
' 45 While direct solution to [3.1] would give the precise 

assignment, the solution of k-dimensional equations directly 

by the set of observations {o-Ji=l M-l}. is not valid for ^ k is 100 com P lex ,ime churning for P?***- 

Because the matrix is sparse the list of cost elements is Thus " * c P resent lnventlon 11115 P roblem 

stored as a packed list and then for each dimension of the 50 The following is a short description of many aspects of the 

matrix, a vector of a variable length list of pointers to the present invention and includes some steps according to the 

cost elements is generated and stored. This organization Pnor art. The first step in solving the problem indirectly is 

means that for a particular observation o. the j* list in the i* to reduce complexity of the problem by the previously 

vector will be a list of pointers to all hypotheses in which o « known ^cussed Ugrangian Relaxation technique. 

- y participates. This structure is further explained in the According to the Lagrangian Relaxation technique, the 

* f. • . v . . absolute requirement that each point is assigned to one and 

following section dealing with problem partitioning. t , • , j . c 6 

& * r tr- * on j v onc 1$ ^1^^ suc h mat f or some one image. 

The objective of the assignment solving process is to points can be assigned to more than one track. However, a 

select from the set of all possible combinations of track ^ penalty based on a respective Lagrangian Coefficient u* is 

extensions a subset that satisfies two criteria. First each added to the cost function when a point in the image is 

point in the subset of combinations should be assigned to assigned to more than one track. The Lagrangian Relaxation 

one and only one track and therefore, included in one and technique reduces the complexity or "dimension" of the 

only one combination of the subset, and second, the total of formulation of the assignment problem because constraints 

the scoring sums far the combinations of the subset should 65 on one observation set are relaxed. Thus, the Lagrangian 

be minimiz ed. This yields the following M-dimensional Relaxation is performed iteratively to repeatedly reduce the 

equations where k=M: dimension until a 2-dimensional penalized cost function 
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problem results as in Problem Formulation |2.1|. This 

2- dimensional problem is solved Chen directly by a previ- 
ously known technique such as Reverse Auction. The penal- 
ized cost function for the 2-dimensional problem defines a 
valley or convex shaped surface which is a function of 
various sets of {u*lk=3, ... ^1} penalty values and one set 
of assignments for the points in two dimensions. That is. for 
each particular u 3 there is a corresponding 2-dimensional 
penalized cost function problem and its solution. Note that 
the solution of the 2-dimensional penalized cost function 
problem identifies the set of assignments for the particular u 3 
values that nunimize the penalized cost function. However, 
these assignments are not likely to be optimum for any 
higher dimensional problem because they were based on an 
initial arbitrary set of u k values. Therefore, the next step is 
to determine the optimum assignments for the related 

3- dimensional penalized cost function problem. There exists 
a 2-dimensional hill shaped function <t> 



10 



15 



18 



function 'I' 



gradient of the t 



results in a new 3-dimensional problem and requires the 
2-dimensional hill climbing based upon new 2-dimensional 
problems. At the peak of the 3-dimensional *F to a complete 
solution. This process is repeated iteratively until the u* 
values that result in a corresponding solution at the peak of 
20 the highest order <t> 



u k \k>3) 



25 



values originally assigned. Then, the gradient of the hill- 
shaped 0 



30 



35 



previously selected for the one point on the hill 
(corresponding to the minimum of the penalized cost <I> 



U 3 



40 



which the corresponding problem will result in the peak of 
the 



45 



characteristic. The next task is to recover the solution of the 
proper u 4 values for the 4-dimensional problem. The fore- 
going hill climbing process will not work again because the 
foregoing hill climbing process when required to locate the 
4-dimensional u 4 values for the peak of <t> 3 exact definition 
of the 4> 3 case of O 2 4> 3 

the iteration can result in a greater than approximation of <I> 3 



Coefficient u M penalty values initially equal to zero. The 
Lagrangian Coefficients associated with the M' h constraint 



50 



set are a N A 



55 



element vector. The reduction of this 



M -dimensional problem in step 322 to a M-l dimensional 
problem uses the two step process described above. First, a 
penalty based on the value of the respective u M coefficient 
is added to the cost function when a point is assigned to more 
than one track and then the resultant cost function is mini- 
mized. However, during the first iteration, the penalty is zero 
because all u M values are initially set to zero. Second, the 
requirement that no point from any image can be assigned to 
more than one track is relaxed for one of the images. In the 
60 extreme this would allow a point from the relaxed image to 
associate with every track. However, the effect of the 
previous penalty would probably mean that such an asso- 
ciation would not minimize the cost The effect of the two 
steps in combination is to remove a hard constraint while 
65 adding the penalty to the cost function so that it operates like 
a soft constraint. For step 322 this two step process results 
in the following penalized cost function problem with k=M: 
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{a) s Minimize: 4* x u 
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•continued 



'k'° L'i* 0 



(£>) Subje 



ject to: \ .■■X z fi "U = 1 il = L - N| 



for ij=U...Nj 
and j = 2, ... k - 1 

«f,.„/ 4 etO. 1) for all 



*fj.„* 4 = Min; arg minimize]^ tii - 
for (i,, 



Jt-i 



0,1 N k \ 

...,ii-])*(0, ...,0) 



[33] 



: £min|0,4.. i0 -«f t | 



Subject to: = 1 



1 2 T,z«,-> 



-0 .^1.0/^,-0 



10 



15 



for 1, ...Nj 
and y = 2, . . . A - 1 



ij*0 - £ - 



1 ii-i = 1 Ai-! 

€(0.1| for all it... ijt-i 
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Because the constraint on single assignment of elements 
from the last image has been eliminated, a M-l dimensional 
problem can be developed by eliminating some of the 
possible assignments. As shown in Equations [3.3]. this is 
done by selecting the smallest cost clement from each of the 
M rh axis vectors of the cost matrix. Reduction in this manner 
yields a new. lower order penalized cost function defined by 
Equations [3.3] which has the same minimum cost as does 
the objective function defined by [3.2] (a) above. 

The cost vectors are selected as follows. Define an index 
array m* f . . . t k j and new cost array c*" 1 ,- • ■ i k _ x by: 



20 Assignment Problem [3.1] and Problem Formulation [3.4] 
differ only in the dimension M vs. M-l. respectively. An 
optimal solution to [3.4] is also an optimal solution to 
equation [3.2]. This relationship is the basis for an iterative 
sequence of reductions indicated by steps 320-332 through 
330-332 and 200-204 in which the penalized cost function 
problem is reduced to a two dimensional problem As these 
formula will be used to describe the processing at all levels, 
the lowercase k is used except where specific reference to 
30 the top level is needed. In step 206. the 2-dimensional 
penalized cost function is solved directly by the prior art 
Reverse Auction technique. Each execution of 206 produces 
two results, the set of z z values that minimize the problem 
and the cost that results v 2 when these z values are substi- 
35 tuted into the objective function [3.4] (a). 

In step 208. according to the prior art. solution z 2 values 
are substituted into the 2-dimensional derivative of the <J> 2 
U should be adjusted so as to perform the hill climbing 
40 function. As was previously described the objective is to 
produce a sequence of u 3 t values which ends when the U 3 p 



value in the domain of the peak of the <I> 



45 



50 



The resulting M-l dimensional problem is (where te=M): 

[3.4] 



55 



to the peak of <P 2 flow moves to step 214 which will attempt 
to recover the 3 -dimensional solution as previously 
described. When further adjustment is required then the flow 
progresses to step 212 and the new values of u* are com- 
puted. At the 2-dimensional level the method of the prior art 
could be used for the hill climbing procedure. However, it is 
not practical to use this prior art hill climbing technique to 
determine the updated Lagrangian Coefficients u* or the 
Max on the next (or any) higher order <J> 
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it is not possible to use the prior art hill climbing to hill climb 
to the peak of <t> k definition of <Z>* functions. As the solution 
of <t> A i in that higher order u values are not yet optimized, its 
value can be larger than the true peak of <t>k lower bound on 
5 v A and it can not be used to recover the required solution. 



m 
o 



10 



values which result in z values that are closer to the proper 
z values for the highest k-dimension. This hill climbing 
process (which is different than that used in the prior art of 
FIG. 4 for recovering only the three dimensional solution) is 
used iteratively at all lower dimensions k of the 
M-dimensional problem (including the 3 -dimensional level 
where it replaces prior art) even when k is much larger than 
three. FIG. 6 helps to explain this new hill climbing tech- 
nique and illustrates the original k-dimensional cost function 
v k of Problem Formulation |3.1|. However, the actual 
k-dimensional cost surface v* defined by 1 3.1 1 (a) comprises 
scalar values at each point described by k-dimensional 
vectors and as such can not be drawn. Nevertheless, for 
illustration purposes only. FIG. 6 ignores the reality of 
ordering vectors and illustrates a concave function v*(z*) to 
represent Equation 3.1. The v* surface is illustrated as being 
smooth to simplify the explanation although actually it can 
be imagined to be terraced. The goal of the assignment 
problem is to find the values of z*; these values minimize the 
k-dimensional cost function v*. 

For purposes of explanation, assume that in FIG. 6, k=4 
(the procedure is used for all k^3). This problem is reduced 
by two iterations of Lagrangian Relaxation to a 
2-dimensional penalized cost function <s>J (k ~ l)i - This cost 
function, and all other cost functions described below are 
also non-smooth and continuous but are illustrated in FIG. 
6 as smooth for explanation purposes. Solving the tyf k ~ l)t 
problem results in one set of z 2 assignments and the value of 
fyk-i) at the point u 2 ,-,. A series of functions tyj k ~ l \ . . . . 
V each generated from a different u 3 are shown. The 
particular series illustrates the process of locating the peak 
of <tV*-iv The 2-dimensional penalized cost functions 
<t>/*~ . • . . « <t>i ( * -1> ' can be solved directly. Each such 
solution provides the information required to calculate the 
next u 3 value. Each iteration of the hill climbing improves 
the selection of u 3 values, i.e., yields <|> 2 solution is closer to 
those at the solution of the 4> 3 The result of solving 
is values that are on both 



and <t> A <t> k 
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Consequently, the z* values at the peak of <P k greatest lower 
bound on the cost of the relaxed problem), can be substituted 
into the k-dimensional penalized cost function to determine 
the proper assignments. Consequently, the present invention 
attempts to find the maximum on the <t> k surface. However. 



15 



and [3.4] were solved at all possible values of u* the function 

than the v k surface except at the peak as described in 
Equation 3.5 and its maximum occurs where the v k surface 
is minimum. Because the <Mc the <t> k <t> k v k surface. The <t> k 
minimization problem described by |3.1| (a). Let z* be the 
unknown solution to Problem Formulation (3.1] and note 
that: 



Instead, the present invention defines all auxiliary func- 
tion 4\ dimensional penalized cost function problem, lower 
order z* values and u* values determined previously by the 
reduction process. The H t k <t> k and its gradient is used for hill 
climbing to its peak. The z* values at the peak of the H' k into 
Problem Formulation |3.2] to determine the proper assign- 
ments. To define the H* k explicitly makes the function <t> k U* 
order sets of Lagrangian Coefficients with the expanded 
notation: <b k U* U* +1 U* 4* is defined recursively, using the 



20 *j(« 3 )=h + X«?j =«i>}(« 3 ;« 4 «*> 

'3=0 



|3.6|<a> 



25 where v 2 is the solution value for the most recent 
2-dimensional penalized cost function problem. For k>3 



I3.6|(b) 
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From the definition of <I> A k compared with Problem Formu- 
lation I3.2|): 



40 



45 it follows that: 



i 3 -0 
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and with that Equation 3.5 is extended to: 



55 



4\(u J , . . . , u^ l ;u*)§O t (u*;u wl u M Sv A (u*)Sv i (u*) [3.71 



This relationship means that either <& k V k hill climbing to 
60 update the Lagrangian Coefficients u. <t> k the preferred 
choice, however it is only available when the solution to 
Problem Formulation 1 3-21 is a feasible solution to Problem 
Formulation 13.1) (as in hill climbing from the second to 
third dimension which is why prior art worked). For sim- 
65 plicity in implementation, the function M\ that it equals <t> k 
therefore always used, even for hill climbing from the 
second dimension. The use of the 4* 
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climbing/peak detection are to determine the gradient of the 
*¥ k and then move up the *¥ k increasing portion of the 
gradient. As shown in FIG. 6. any upward step on the % V k the 
4^- U* a a is closer to the ideal set of u* values which 
correspond to the minimum of the v k function. While it is 
possible to iteratively step up this 4\ and then determine if 
the peak has been reached, the present invention optimizes 
this process by determining the single step size from the 
; point at the ininimum of $ k that will jump to the 
: and then calculating the u* values at the peak. Once the 
r i. "u" suat the peak at x ¥ k determined, then the u* values can 
Sstituted into Problem Formulation [3.2] to determine 
the proper assignment. (However, in a more realistic 
example, where k is much greater than three, then the u* 
values at the peak of the 4* along with the lower order u* 
values and those assigned and yielded by the reduction steps 
are used to define the next higher level 



determine the gradient of each 4* 
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The modification to Wolfe* s algorithm uses the information 
generated for % ¥ k U* 3 U**" 2 U"" 1 for calculating the sub- 
gradients. The definition of a subgradient v of 4^ U* 3 U** 
subdifferential set defined as: 

5 54' 4 (u>=}v€R^llOJ' t (U* 3 U^IT)-^" • • • . 

U*- 1 ;U*))Sv T CLT-U*) WeR *1} 

(where v r is the transpose of v) Next, a subgradient vector 
is determined from this function. If z* is the solution of 
10 Problem Formulation |3.2], then differentiating 4\ U ? U*" 1 
U* U* a evaluating the result with respect to the current 
selection matrix z* yields a subgradient vector: 
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"bundle") are deterrnined at several points on the 4* t in the 
region of the starting point (Le., minimum of <t>* ) and then 
averaged. Statistically, the averaged result should point 
toward the peak of the % ¥ k Subgradient Algorithm (Wolfe75, 
Wolfe79) for nuriimization was previously known in another 
environment to determine a gradient of a non-smooth sur- 
face using multiple subgradients and can be used with 
modification in the present invention. Wolfe's algorithm is 
further described in "A Method of Conjugate Subgradients 
for Minirnizing Nondifferentiable Functions" page 147-173 
published by Mathematical Programming Study 3 in 1975 
and "Finding the Nearest Point in a Polytope" page 128-149 
published by Mathematical Prograrnming Study 11 in 1976. 



The iterative nature of the solution process at each dimen- 
sion yields a set of such subgradients. Except for a situation 
described below where the resultant averaged gradient does 
not point toward the peak, the most recent set of such 
subgradients are saved and used as the "bundle" for the peak 
finding process for this dimension. For example, at the k 
level there is a bundle of subgradients of the *¥ k the mini- 
mum of the <)>*. surface determined as a result of solving 
Problem Formulation [3.1] at all lower levels. This bundle 
can be averaged to approximate the gradient. Alternately, the 
previous bundle can be discarded so as to use the new value 
to initiate a new bundle. This choice provides a way to adjust 
the process to differing classes of problems, i.e.. when data 
is being derived from two sensors and the relaxation pro- 
ceeds from data derived from one sensor to the other then the 
prior relaxation data for the first sensor could be detrimental 
to performance on the second sensors data. 
L2.3. Step Size and Termination Criterion 
After the average gradient of the 4* k determined, the next 
step is to determine a single step that will jump to the peak 
of the 4 / Jt is to first specify an arbitrary step size, and then 
calculate the value of 4\ gradient. If the *¥ k this probably 
means that the step has not yet reached the peak. 
Consequently, the step size is doubled and a new x ¥ k value 
determined and compared to the previous one. This process 
is repeated until the new 4** previous one. At that time, the 
step has gone too far and crossed over the peak. 
Consequently, the last doubling is rolled-back and that step 
size is used for the iteration. If the initial estimated 4\ step 
size is decreased by 50%. If this still results in a smaller 4\ 
previous step size is decreased by 25%. The following is a 
more detailed description of this process. 

With a suitable bundle of subgradients determined as just 
described. Wolfe's algorithm can be used to determine the 
effective subgradient (d) and the Upgraded value U*^. 
From the previous iteration, or from an initial condition, 
there exists a step length value (t). The value. 

is calculated as an estimate of u*^. To determine if the 
60 current step size is valid the we evaluate 4* fc u 3 u*~ 2 u + If the 
result represents an improvement then we double the step 
size. Otherwise we halve the step size. In either case a new 
u + is calculated. The doubling or halving continues until the 
step becomes too large to improve the result or until it 
65 becomes small enough to not degrade the result The result- 
ing suitable step size is saved with d as part of the subgra- 
dient bundle. The last acceptable u + is assigned to u*^. 
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for / = 0 N k 

and y = 0,...,Afo 
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If the resulting 2-dimensional assignment problem, 

% % 

Minimize: ^ ^ A ^p w p 
/-o /*o 

Subject to: J = 1 Wo 

** 

M Wo 

/«o 

w M e (0, 1| J = 0. ... , / = 0, ... , N k 
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Three distinct criterion are used to determine when u k j is 
close enough to u*: 

1. The Wolfe's algorithm criterion of d=0 given that the 
test has been repeated with the bundle containing only 
the most recent subgradient. 5 

2. The difference between the lower bound <t> k u k upper 
bound v k (z*, u*k) being less than a preset relative 
threshold. (Six percent was found to be an effective 
threshold for radar tracking problems.) 

3. An iteration count being exceeded. 
The use of limits on the iteration are particularly effective 

for iterations at the level 3 through n-1 levels, as these 
iterations will be repeated so as to resolve higher order 
coefficient sets. With limited iterations the process is in 
general robust enough to improve the estimate of upper 15 
order Lagrangian Coefficients. By limiting the iteration 
counts then the total processing time for the algorithm 
becomes deterministic. That characteristic means the pro- 
cess can be effective for real time problems such as radar, 
where the results of the last scan of the environment must be 20 
processed prior to the next scan being received. 
L2.4. Solution Recovery 

The following process determines if the u* values at what 
is believed to be the peak of the H* k approximate the u* 
values at the peak of the corresponding <t> k function. This is 25 
done by determining if the corresponding penalized cost 
function is minimized. Thus, the u* values at the peak of the 
% V k Formulation |3.2] to determine a set of z assignments for 
k-1 points. During the foregoing reduction process. Problem 
Formulation |3.4| yielded a tentative list of k selected z 30 
points that can be described by their indices as: {(F,. . . 
F Jk _ 1 )} A ^_ 1 , where N 0 is the number of cost elements 
selected into the solution. One possible solution of Problem 
Formulation [3.1] is the solution of Problem Formulation 
[3.2] which is described as {(rV . . r^m^. . . i*.,)}" 0 ^ 35 
with m* . . . i A _! as it was defined in Problem Formulation 
[3.3]. If this solution satisfies the k-th constraint set then it 
is the optimal solution. 

However, if the solution for Problem Formulation [3.2 1 is 
not feasible (decision 355), then the following adjustment 40 
process determines if a solution exists which satisfies the 
k-th constraint while retaining the assignments made in 
solving Problem Formulation I3.4J. To do this, a 
2-dimensional cost matrix is defined based upon all obser- 
vations from the k-th set which could be used to extend the 45 
relaxed solution. 
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has a feasible solution then the indices of that solution map 
to the solution of Problem Formulation [3.1] for the 



k-dimensional problem. The first index in each resultant 
solution entry is the pointer back to an element of the \(¥ r 
. . if k _ l m k i . . . i*_i)r >=i list. That element supplies the first 
k-1 indices of the solution. The second index of the solution 
to the recovery problem is the k"' index of the solution. 
Together these indexes specify the values of z* that solve 
Problem Formulation |3.1| at the k" 1 level. 

If Problem Formulation |3.9| does not have a feasible 
solution then the value of uki which was thought to represent 
u*^ is not representative of the actual peak and further 
iteration at the k th level is required. This decision represent 
the branch path from step 214 and equivalent steps. 

1.3. Partitioning 

A partitioning process is used to divide the cost matrix 
that results from the Scoring step 154 into as many inde- 
pendent problems as possible prior to beginning the relax- 
ation solution. The partitioning process is included with the 
Problem Formulation Step 310. The result of partitioning is 
a set of problems to be solved, i.e., there will be p, problems 
that consist of a single hypothesis, p 2 problems that consist 
of two hypothesis, etc. Each such problem is a group in that 
one or more observations or' tracks are shared between 
members of the group. 

In partitioning to groups no consideration is given to the 
actual cost values. The analysis depends strictly on the basis 
of two or more cost elements sharing the same specific axis 
of the cost matrix. In a 2-dimensional case two cost elements 
must be in the same group if they share a row or a column. 
If the two elements are in the same row. then each other cost 
clement that is also in the same row, as well as any cost 
elements that are in columns occupied by members of the 
row must be included in the group. The process continues 
recursively. In literature it is referred to as "Constructing a 
Spanning Forest." The k-dimensional case is analogous to 
the 2-dimensional case. The specific method we have incor- 
porated is a depth first search, presented by Aho, Hopecraft, 
and Ullman, in "Design and Analysis of Computer 
Algorithms," section 5.2, published by Addison-Westley, 
Mass. 

The result of partitioning at level M is the set of problems 

described as \P 0 \i=l p y and j=l N}, where N is 

the total number of hypothesis. The problems labeled {P a li= 

1 pj } are the cases where their is only one choice for 

the next observation at each scan and that observation could 
be used for no other track, i.e.. it is a single isolated track. 
The hypothesis must be in included in the solution set and 
no further processing is required. 

As hypothesis are constructed the first clement is used to 
refer to the track id. Any problem in the partitioned set which 
does not have shared costs in the first scan represent a case 
where a track could be extended in several ways but none of 
the ways share an observation with any other track. The 
required solution hypothesis for this case is the particular 
hypothesis with the maximum likelihood. For this case all 
likelihoods are determined as was described in Scoring and 
the maximum is selected. 

In addition to partitioning at the M level, partitioning is 

applied to each subsequent level M-l 2. For each 

problem that was not eliminated by either by the prior 
selection partitioning is repeated, ignoring observations that 
are shared in the last set of observations. Partitioning rec- 
ognizes that relaxation will eliminate the last constraint set 
and thus partitioning is feasible for the lower level problems 
that will result from relaxation. This process is repeated for 
all levels down to k=3. The full set of partitionings can be 
performed in the Problem Formulation Step 310, prior to 
initiating the actual relaxation steps. The actual 
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2-dimensional solver used in step 206 includes an equivalent 
process so no advantage would be gained by partitioning at 
the te=2 level. 

There are two possible solution methods for the remaining 
problems. "Branch and Bound" as was previously described, 
or the relaxation method that this invention describes. If any 
of the partitions have 5-10 possible tracks and less than 50 
to 20 hypotheses, then the prior art 4t Branch and Bound" 
algorithm generally executes faster than does the relaxation 
due to its reduced level of startup overhead. The "Branch 
and Bound" algorithm is executed against all re m a inin g M 
level problems that satisfy the size constraint For the 
remainder the Relaxation algorithm is used. The scheduling 
done in Problem Formulation allows each Problem Formu- 
lation [3.2 1 cost matrix resulting from the first step of 
relaxation to be partitioned. The resulting partitions can be 
solved by any of the four methods, isolated track direct 
inclusion, isolated track tree evaluation, small group 
"Branch and Bound" or an additional stage of relaxation as 
has been fully described. 

The partitioning after each Level of Lagrangian Relax- 
ation is effective because when a problem is relaxed, the 
constraint that requires each point to be assigned to only one 
track is eliminated (for one image at a time) . Therefore, two 
tracks previously linked by contending for one point will be 
unlinked by the relaxation which permits the point to be 
assigned to both tracks. The fruitfulness of partitioning 
increases for lower and lower levels of relaxation. 

The following is a more detailed description of the 
partitioning method. Its application at all but the M level 
depends upon the relaxation process described in this inven- 
tion. The recursive partitioning is therefore a distinct part of 
this invention. The advantage of this method is greatly 
enhanced by the sparse cost matrix resulting from tracking 
problems. However the sparse nature of the problem 
requires special storage and search techniques. 

A hypothesis is the combination of the cost element c„ the 
selection variable z n and all observations that made up the 
potential track extension. It can be written as, h n ={c n . z„. 

{o n =okiik=l Mie{ 1 Nk}}}. i.e.. cost, selection 

variable and observations in the hypothesis. ne{ 1 N} 

where N is the total number of hypotheses in the problem. 
While the cost and assignment matrices were previously 
referenced, these matrices are not effective storage mecha- 
nisms for tracking applications. Instead the list of all hypoth- 
esis and sets of lists for each dimension that reference the 
hypothesis set are stored. The hypothetical set in list form is: 

For each dimension k=l . . . M there exists a set of lists, with 
each list element being a pointer to a particular hypothesis: 

where is number of hypothesis containing the i-th 
observation from scan k. This structure is a multiply linked 
list in that any observation is associated with a set of pointer 
to all hypothesis it participates in. and any hypothesis has a 
set of pointers to all observations that formed it. (These 
pointers can be implemented as true pointers or indices 
depending upon the particular programming language 
utilized.) 

Given this description of the matrix storage technique 
men the partitioning technique is as follows: Mark the first 
list Ljt ; and follow out all pointers in that list to the indicated 

hypothesis h^ for i=l Nki. Mark all located 

hypothesis, and for each follow pointers back the particular 
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L 0k for k=l M. Those L's if not previously marked get 

marked and also followed out to hypothesis elements and 
back to L's. When all such L's or h's being located are 
marked, then an isolated sub-problem has been identified. 
5 All marked elements can be removed from the lists and 
stored as a unique problem to be solved. The partitioning 
problem then continues by again starting at the first residual 
L set. When none remain, the original problem is completely 
partitioned. 

Isolated problems can result from one source track having 
multiple possible extensions or from a set of source tracks 
contending for some observations. Because one of the 
indices of k (in our implementation it is k=l) indicate the 

15 source track then it is possible to categorize the two problem 
types by observing if the isolated problem includes more 
than one L-list from the index level associated with tracks. 
L4. Subsequent Track Prediction 
As noted above, after the points are assigned to the 

20 respective tracks, some action is taken such as controlling 
take-offs and landings to avoid collision, advising airborne 
aircraft to change course, warning airborne aircraft of an 
adjacent aircraft in a commercial environment, or aiming 

25 and firing an anti-aircraft (or anti-missile) missile, rocket or 
projectile, or taking evasive action in a military environ- 
ment Also, the tracking can be used to position a robot to 
work on an object. For some or all of these applications, 
usually the tracks which have just been identified are 

30 extended or extrapolated to predict a subsequent position of 
the aircraft or other object. The extrapolation can be done in 
a variety of prior art ways; for example, straight line 
extensions of the most likely track extension hypothesis, 
parametric quadratic extension of the x versus time, y versus 

35 time, etc. functions that are the result of the filtering process 
described earlier, least square path fitting or Kalman filtering 
of just the selected hypothesis. 
The process of gating, filtering, and gate generation as 

^ they were previously described require that a curve be fit 
through observations such that fit of observations to the 
likely path can be scored and further so that the hypothetical 
target future location can be calculated for the next stage of 
gating. In the implementation existing, quadratic functions 

45 have been fit through the measurements that occur within the 
window. A set of quadratics is used, one per sensor mea- 
surement. To calculate intercept locations these quadratics 
can be converted directly into path functions. Intercept times 
are then calculated by prior methods based upon simulta- 

50 neous solution of path equations. 

The use of the fitted quadratics is not as precise as more 
conventional filters like the Kalman Filter. They are however 
much faster and sufficiently accurate for the relative scores 
required for the Assignment problem. When better location 

55 prediction is required then the assignment problem is 
executed to select the solution hypothesis and based upon 
the observations in those hypothesis the more extensive 
Kalman filter is executed. The result is tremendous compu- 

6q tation savings when compared with the Kalman Filter being 
run on all hypothetical tracks. 

Based on the foregoing, apparatus and methods have been 
disclosed for tracking objects. However, numerous modifi- 
cations and substitutions can be made without deviating 

65 from the scope of the present invention. For example, the 
foregoing 4* u* values are used to eliminate the higher than 
approximation characteristic of the <t> k 
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be partitioned into tracks and false alarms. In the formula- 

* tion of track extensions a moving window over time of 

observations sets is used. The observation sets will be 

measurements which are: (a) assigned to existing tracks, (b) 

<P t 5 designated as false measurements, or (c) used for initiating 

tracks. However, note that alternative data objects instead of 

observation sets may also be fused such as in sensor level 
be as efficient it seems likely that the method would con- whefein each sensof fonns Uacks from its own 

verge. Therefore the invention .has i been disclosed by way of measurements and then me tracks from the sensors ^ fused 
example and not limitation, and reference should be made to io ^ ^ as ^ ^ ^ ^ wm 

the following claims to determine the scope of the present . ... . 

6 r appreciate, both embodiments of the present invention may 

invention. ^ * * L . *. f^**- 

be used for this type of data fusion. 

II. AN ALTERNATIVE EMBODIMENT OF THE The data assignment problem considered presently is 
MULTI-DIMENSIONAL ASSIGNMENT 15 represented as a case of set partitioning defined in the 

SOLVING PROCESS following way. First, for notational convenience in repre- 

The following description provides an alternative second senting tracks, a zero index is added to each of the index sets 

embodiment of the multi-dimensional assignment solving in 1 5.2 1. a gap filler Z* 0 is added to each of the data sets Z(k) 

process of section 1.1. However, in order to clearly describe in |5.2|. and a "track of data" is defined as (z 1 ^ z N if ) 

the present alternative embodiment, further discussion is 20 where i* and z\ 4 can now assume the values of 0 and z*^ 

first presented regarding the data assignment problem of respectively. A partition of the data will refer to a collection 

partitioning measurements into tracks and false alarms and, of tracks of data or track extensions wherein each observa- 

in particular, the representation of multi-dimensional assign- tion occurs exactly once in one of the tracks of data and such 

ment problems. that all data is used up. Note that the occurrence of the gap 

ELL Formulation of the Assignment Problem The goal of 25 filler is unrestricted The gap filler z* 0 serves several pur- 

this section is to explain the formulation of the data asso- poses in the representation of missing data, false 

ciation problems, and more particularly multi-dimensional observations, initiating tracks, and terminating tracks. Note 

assignment problems, that govern large classes of problems that a "reference partition" may be defined which is a 
in centralized or hybrid centralized-sensor level multisensor/ ^ partition in which all observations are declared to be false, 

multitarget tracking The presentation is brief; technical Next undef a iate indepen dence assumptions it can 

detads are presented for both track initiation and mamte- be shown that . 
nance in A. B. Poore, Multi-dimensional assignment formu- 
lation of data association problems arising from multitarget 

tracking and multisensor data fusion, Computational Opti- P\X = y\2r) 2 . H pi ^ / ^ 

mizationandApplications,3(1994).pp.27-57,fornonma- f\Y=yfi\z*) „ l „ , ^ )cy 
neuvering targets and in Ingemar J. Cox, Pierre Hansen, and 
Beta Julesz, eds.. Partitioning Data Sets, DIMACS Series in 

Discrete Mathematics and Theoretical Computer Science. where L*^ i M is a likelihood ratio containing probabilities 

American Mathematical Society. Providence, R.L, v. Feb. for detection, maneuvers, and termination as well as prob- 

19, 1995, pp. 169-198 for maneuvering targets. Note that the ability density functions for measurement errors, track ini- 

present formulation can also be modified to include target tiation and termination. Then if c*, j^-lnL*^ iM , 

features (e.g., size and type) into the scoring process 154. 1 ' " 

The data assignment problems for multisensor and mul- i j\ Y \z M ) I ^ I 5 - 4 ! 

titarget tracking considered here are generally posed as that 45 ^ fZ « I = X c 'j 

of maximizing the posterior probability of the surveillance t'i.»Jn>«r 
region (given the data) according to 



w . . [f(r = y|2 M ) 
Maximize 



50 ^ 



Using [5.31 and the zero-one variable z\ t^l if (i A 

) 

€ y I5.5| 



where Z M represents M data sets, y of the data (and thus 
induces a partition of the data), V* the finite collection of all 
such partitions, T 55 

P(rYlZ M ) is the posterior probability of a partition y true y ,j> y _ Vc, ( 

given the data Z M . The term partition is defined below. /|-o i^-o 1 1 

Consider M observation sets Z(k) (k=l N) each of M 

N* observations {z*^}^-!, and let Z M denote the cumulative 60 y V . f \ u iy = i, 

data set defined by £t 1 " 

(o X — X X — X *!•«'* = l > 

respectively. In multisensor data fusion and multitarget fa i L ^uti' 0 
tracking the data sets Z(k) may represent different classes of 65 

objects, and each data set can arise from different sensors. for '* = 1 Afc and * « 2. .... J# - 1. 

For track initiation the objects are measurements that must 
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-continued 



(d) 
(e) 



-1, '« = 1 

for all i| i'm 



where Cq 0 is arbitrarily defined to be zero. Here, each 
group of sums in the constraints |5.5] (b)-[5.5] (e) repre- 
sents the fact that each non-gap filler observation occurs 
exactly once in a t4 track of data." One can modify this 
formulation to include multiassignments of one, some, or all 
the actual observations. The assignment problem |5.5] is 
changed accordingly. For example, if z k ik is to be assigned no 
more than, exactly, or no less than n* t times, then the "=1** 
"in the appropriate one of the constraints |5.5 | (b)— [5.5] (d) 
is changed to =. ^n*^." respectively. 

Expressions for the likelihood ratios L* j - . . i M such as [5.5] 
are well known. In particular, discussions of these ratios may 
be found in: A.B. Poore. Multi-dimensional assignment 
formulation of data association problems arising from mul- 
titarget tracking and multisensor data fusion. Computational 
Optimization and Applications. 3 (1994). pp. 27-57; and 
Ingemar J. Cox, Pierre Hansen, and Beta Julesz. eds.. 
Partitioning Data Sets. DIMACS Series in Discrete Math- 
ematics and Theoretical Computer Science. American Math- 
ematical Society. Providence. R.L. v. 19, 1995. pp. 169-198. 
Furthermore, the likelihood ratios are easily modified to 
include target features and to account for different sensor 
types. Also note that in track initiation, the M observation 
sets provide observations from M sensors, possibly all the 
same. Additionally note that for track maintenance, a sliding 
window of M observation sets and one data set containing 
established tracks may be used. In this latter case, the 
formulation is the same as above except that the dimension 
of the assignment problem is now M+l. 

H.2. Overview of the New Lagrangian Relaxation Algo- 
rithms 

Having formulated an M-dimensional assignment prob- 
lem [5.5]. we now turn to a description of the Lagrangian 
relaxation algorithm. Subsection H.2.1 below presents many 
of the relaxation properties associated with the relaxation of 
an if dimensional assignment problem to an m-dimensional 
one via a Lagrangian relaxation of n-m sets of constraints 
wherein m<n<M and preferably in the present invention 
embodiment n-m>l. Although any n-ro constraint sets can 
be relaxed, the description here is based on relaxing the last 
n-m sets of constraints and keeping the first m sets. Given 
either an optimal or suboptimal solution of the relaxed 
problem, a technique for recovering a feasible solution of the 
n-dimensional problem is presented in subsection IL2.2 
below and an overview of the Lagrangian relaxation algo- 
rithm is given in subsection H.2.3. 

The following notation will be used throughout the 
remainder of the work. Let M be an integer such that M^3 
and let ne{3, . . . M}. The n-dimensional assignment 
problem is 



*1 #m 

(a) Minimize: v m {z)^ £ ... £<3i„ .^...4, 
1*1 "m 

<b) Subjectto: £ ... B '1 0 1 



[6.1) 



(c) 
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(e) 
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-continued 

Z- Z Z~Zw- 

for k = 1 N k 

and Jfc = Z ... , n - 1, 

Z- Z N - 



eta Ik 



for all i\ i„ 



15 To ensure that a feasible solution of [6.1] always exists, all 
variables with exactly one nonzero index (i.e.. variables of 
the form 2T 0 0 , t0 0 for i^^O) are assumed free to be 
assigned and the corresponding cost coefficients are well- 
defined. 

20 

II 2.1. The Lagrangian Relaxed Assignment Problem 

The n-dimensional assignment problem [6.1] has n sets of 
constraints. A (N*+l)-dimensional multiplier vector associ- 

25 ated with the k-th constraint set will be denoted by u*=(u* 0 . 

u*, u k M ) T with u* o =0 and k=l n. The 

n-dimensional assignment problem [6.1] is relaxed to an 
m-dimensional assignment problem by incorporating n-m of 
the n sets of constraints into the objective function |6.1] (a). 

30 Although any constraint set can be relaxed, the description 
of the relaxation procedure for [6.1] will be based on the 
relaxation of the last n-m sets of constraints. The relaxed 
problem is 
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: Minimize 



|6.2] 



•4,-1 



,£...2k..., + z<k....-zz< 



Subject to: 



♦2-0 



< 1 Afi. 



for u = 1 N k and k « 2, .... m, 

^...^.elO, 1} for all i, i„. 

60 wherein we have added the multiplier u* o =0 for notational 
convenience. Thus, the multiplier R with u* o =0 is fixed. 

An optimal (or suboptimal) solution of [6.2] can be 
constructed from that of an m-dimensional assignment prob- 
65 lem. To show this, def ine f or each (ij i„) an index 

Gm+r • • • J«) ^Orn^lOl im) i J) Httd Sk 

new cost function c" 1 ^. • • ' m by 
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-continued 

OU+i ji,)=argmin I 6 - 3 ! h^..^ =0 if ... O * U«fi» • ■• Ju) and 



k = 0, . , ty. and Jt = m + 1 n 



II 

Z "*A for( " i-)*(0,..-.O) 



(ii i»)*(0.....0) 

■ 



10 



< *-° aB S -■- 4 ^ Mi ^ lim { 0 ^ c S»-«-*i--.4- + Z "5*j 

'* + i=o i„-o * , " B+1 ^hen w « i s a feasible solution of the Lagrangian relaxed 

problem |6.2] and 

(If (j^i j„n) is not unique, choose the smallest such 15 

j^j, amongst those (n-m)-tuples with the same j m+1 choose <u, f m + x n ^ w«;ir +, t . . . t uV^^.i^-ou*^ 

the smallest j^, etc., so that G m+ i j«) * s uniquely . « . . 

defined.) Usi£fte cost coeffid«£ defined in this way. the * « n addl ^ n - j" ,s °P tunal for ' 6 - 4 '- then «" 15 an °P Umal 

following m-dimensional assignment problem is obtained: so u ion o | . f an 

20 o^y- 1 " <i>™» »n-i\^ii"Vw\ 

With the exception of one equality being converted to an 
A *«> inequality, the following Fact is a converse of Fact A.2 and 

s 2_j "' Zj ,fl,z *I has also been shown to be true. 
'i tt0 ' m=0 25 pact Let w" be a feasible solution to problem 16.2] 

and define W m by 



Subject to: > zJJ .../„ = 1, N|, 



'i-o i 4 .i-oi UI -o Hff„. 0 =0 if (i, i w ) = (0,...,0) and 

for = l and Jt = 2, .... * m 



(6.6] 



,» for Hi i«)*(0,...,0) 
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> - Z 3 • '» = 1 

^ h£...o = 1 if <<i U) - (0, .... 0) and 

$ ...i„e{0, 1} foraUi! i„. A 

^-Wfti+t-. + 2-» % ~ 501116 Ib) ' 

40 

As an aside, observe that any feasible solution z" of [6.1] 

yields a feasible solution z m of [6.4] via the construction Then w^isa feasible solution of the problem [6.4] and 

!1 if tf,...^.,^ = 1 for son* k> 
45 If, in addition, w" is optimal for 16.2], then is an optimal 
0 othe,wisc - solution of [6.41, m+1 n 

^(w^;ir +1 vn-ir^l"*^^ and o>_ " 

Thus, the mn-dimensional assignment problem |6.4| has at 

least as many feasible solutions of the constraints as the $Q l * - ■ . «>^w,i£*V>uV 

original problem |6.1 1. n 2 2 The Recovery procedure 

The following Fact A.2 has been shown to be true. It states ^ next o5 j ective is t0 explain a recovery procedure or 
that an optimal solution of Problem Formulation |6.2| can be met hod for recovering a solution to the n-dimensional 
computed from that of Problem Formulation |6.4| pro blem of 1 6.1 1 from a relaxed problem having potentially 
Furthermore, the converse to this fact is provided in Fact 55 substantially fewer dimensions than |6.1|. Note that this 
A.3. Moreover, if the solution of either of these two prob- aspecl of me a | ternauV e embodiment of the present inven- 
lems |6.2| or |6.4| is e-optimal (i.e., the objective function tion is substantially different from the method disclosed in 
associated with the suboptimal solution is within V\ of the the ^ embodiment of the multi-dimensional assignment 
objective function associated with the optimal solution), solving process of section 1.1 of this specification. Further, 
then so is the other. 60 this alternative embodiment provides substantial benefits in 

Fact A.2. Let w" 1 be a feasible solution to problem [6.4] terms of computational efficiency and accuracy over the first 
and define w" by embodiment, as will be discussed. Thus, given a feasible 

(optimal or suboptimal) solution w™ of [6.4] (or w" if 
4, = "ft.. j» if •• ^ = ... Jm) and f 6 5 l Problem Formulation [6.2] is constructed via Fact A.2), an 

65 explanation is provided here regarding how to generate a 

(/| W * (0, .... 0) feasible solution z n of [6. 1] which is close to w" in a sense 

to be specified. We first assume that no variables in [6. 1] are 



5.959.574 



35 



preassigned to zero (this assumption will be removed 
shortly). The difficulty with the solution w" is that it need not 
satisfy the last n-m sets of constraints in [6-1]- Note, 
however, that if is an optimal solution for |6.4] and w" 
(constructed as in Fact A.2) satisfies the relaxed constraints, 
then w n is optimal for 16-1]. The recovery procedure 
described here is designed to preserve the 0-1 character of 
the solution w" of (6.4] as far as possible. That is, if w"V. . . 

i m =l and ilw*0 for at least one 1=1 m. then the 

corresponding feasible solution z* of 1 6.1] is constructed so 
that z n ix . . . i„=l for some (w. . . i„). Note that by mis 
reasoning, variables of the form z n 0 ow . . ^ in can be 
assigned to a value of 1 in the recovery problem only if 
w«q 0 =1. However, variables z n 0 . oc+| . . mim will be treated 
differently in the recovery procedure in that they can be 
assigned 0 or 1 independent of the value w"^ ^ This 
increases the feasible set of the recovery problem, leading to 
a potentially better solution. 

Let {(F A . . . . JfJ}"*^ be an enumeration of indices of 
w" (or the first m indices of wl constructed as in Fact A.2) 

such that w"V . . and & l FJ*fft . . . , 0) . Set 

{f 1 ,OkO 0) for j=0 and define 

-J far M 

if.„4k. + i™* jt = m + l,...,n;; = 0,...,tf 0 . 



Let Y denote the solution of the (n-m+l)-dimensional 
assignment problem: 



[6.8] 



Minimize ^ 

> •■•£>w..<. = 1 ' J- 1 



Subject to 



*0 



for i t = 1, N t 

and k - m + 1 n - 1, 

k - 1 AT., 

y/w*|...t € |0, 1} for an j, 

Hie recovered feasible solution z n of [6.1] corresponding to 
the multiplier set {u m u n } is then defined by 



1 if (ii im) = (f{ for some ^ 

yoO,...,M> and Yji^.i. <= 1 
0 otherwise. 



This recovery procedure is valid as long as all cost coeffi- 
cients c" are defined and all zero-one variables in z" are free 
to be assigned. Note that modifications are necessary for 
sparse problems. If the cost coefficient c% . . l mim¥V . . ^ is well 
defined and the zero-one variable z\. . . Ohh-i- . . i„ is free 
to be assigned to zero or one, then define c"-"** 1 



<* fi{ ■ ? m Wi. • . i„ as in [6.7] with zT*" 1 ^ . . /„ being free 
to be assigned. Otherwise, t**"* 1 ^ . is preassigned to 
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zero. To ensure that a feasible solution exists, we now only 
need ensure that the variables z n ~ m * l J0 0 are free for. 

j^O.1 N 0 . (Recall that all variables of the form 

z n 0 i o are free (to be assigned) and all coefficients of the 

5 form c^ 0 ' v . . 0 are well defined for k=i n.) This 

is accomplished as follows. If the cost coefficient c rt . . . 
F m 0 ... 0 is well defined and z" if j J f m 0 ... 0 is free, then 
define c"-^ 1 * . . . cr*V 2 - ■ - 0> • • 0 with z*-""Vp . . . 0 being 
free. Otherwise, since all variables of the form z* 0 . . (V . . 0 
are known feasible and have well-defined costs, put 



10 



15 



n-m+1 . 

c jO...O 



II.3. The Multi-Dimensional Lagrangian Relaxation 
Algorithm for the n-Dimensional Assignment Problem 

Starting with the M-dimensional assignment problem 
1 6. l.|. i.e.. n=M. the algorithm described below is recursive 

20 in that the M-dimensional assignment problem is relaxed to 
an m-dimensional problem by incorporating (n-m) sets of 
constrains into the objective function using Lagrangian 
relaxation of this set This problem is maximized with 
respect to the Lagrange multipliers, and a good suboptimal 

25 solution to the original problem is recovered using an 
(n-m+1) -dimensional assignment problem Each of these 
two (the m-dimensional and the (n-m+1) -dimensional 
assignment problems) can be solved in a similar manner. 
More precisely, reference is made to FIG. 8 which pre- 

30 sents a flowchart of a procedure embodying the multi- 
dimensional relaxation algorithm, referred to immediately 
above. This procedure, denoted MULH 13 DIM_RELAX in 
FIG. 8. has three formal parameters, n. PROB_ 
FORMULATION and ASSIGNMT_SOLUT10N, which 

35 are used to transfer information between recursive instan- 
tiations of the procedure. In particular, the parameter, n. is an 
input parameter supplying the dimension for the multi- 
dimensional assignment problem (as in [6.1]) which is to be 
solved (i.e.. to obtain an optimal or near-optimal solution). 

40 The parameter. PROB_FORMULATION. is also an input 
parameter supplying the data structure(s) used to represent 
the n-dimensional assignment problem to. be solved. Note 
that PROB_FORMULXTION at least provides access to the 
cost matrix, |c rt ], and the observation sets whose observa- 

45 tions are to be assigned. Additionally, the parameter. 
ASSIGNMT__SOLUTION. is used as an output parameter 
for returning a solution of a lower dimensioned assignment 
problem to an instantiation of MULTI_DIM_RELAX 
which is processing a higher dimensioned assignment prob- 

50 lem. 

A description of FIG. 8 follows. Assuming MULTL_ 
DIM_RELAX is initially invoked with M as the value for 
the parameter n and the PROB_ FORMULATION having a 
data structures) representing an M-dimensional assignment 

55 problem as in [5.51. in step 500 an integer m. 2^m<n is 

chosen such that the constraint sets or dimensions m+1 

n are to be relaxed so that an m-dimensional problem 
formulation as in [6.21 results. In step 504 an initial approxi- 
mation is provided for {u^q u n o56 . Subsequently, in 

60 step 508 the above initial values for {u BH " , 0 ._. • - .u"o} are 

used in an iterative process in determining (u"** 1 u n ) 

which maximizes {^ mn 'f t±l n *eR Ml+l and k=m+l. . . . *} for 
a feasible solution w^ subject to the constraints of Problem 
Formulation [6.2]. Note that by irmimizing O^Y" 1 " 1 " 

65 constraints mat were relaxed are being forced to be satisfied 
and in so doing information is built into a solution of this 
function for solving the input assignment problem in 



5 T 9i 
37 

"PROB_FORMULATION." Further note that a non- 
smooth 10 optimization technique is used here and that a 
preferred method of determining a maximum in step 508 is 
the bundle-trust method of Schramm and Zowe as described 
in H. Schramm and J. Zowe. A version of the bundle idea for 
minimizing a non-smooth function: Conceptual idea, con- 
vergence analysis, numerical results, SIAM Journal on 
Optimization, 2 (1992). pp. 121-152. This method, along 
with various other methods for determining the maximum in 
step 508. are discussed below. 

Also note that for m>2, a solution to the optimization 
problem of step 508 is NP-hard and therefore cannot be 
solved optimally. That is, there is no known computationally 
tractable method for guarantying an optimal solution. Thus, 
there are two possibilities: either (a) allow m to be greater 
than 2 and use auxiliary functions similar to those disclosed 
in the first embodiment of the k-dimensional assignment 
solver 300 in section I to compute a near-optimal solution, 
or (b) make m=2 wherein algorithms such as the forward/ 
reverse auction algorithm of D. P. Bertsekas and D. A. 
Castafion in their paper, A forward/reverse auction algorithm 
for asymrmnetric assignment problems. Computational 
Optimization and Applications. 1 (1992), pp. 277-298, 
provides an optimal solution. 

If option (a) immediately above is chosen, then the 
auxiliary functions are used as approximation functions for 
obtaining at least a near-optimal solution to the optimization 
problem of step 508. Note that the auxiliary functions used 
depend on the value of m. Thus, auxiliary functions used 
when m=3 will likely be different from those for m=4. But 
in any case, the optimization procedure is guided by using 
the merit function, O^y 3 " 2-dimensional assignment prob- 
lem for guiding the maximization process. 

Alternatively, if option (b) above is chosen, then two 
important advantages result: the optimization problem of 
step 508 can be always solved optimally and without using 
auxiliary approximation functions. Thus, better solutions to 
the original M-dimensional assignment problem are likely 
since there is no guarantee that when the non-smooth 
optimization techniques are applied to the auxiliary func- 
tions the techniques will yield an optimal solution to step 
508. Furthermore, it is important to note that without aux- 
iliary functions, the processing in step 508 is both concep- 
tually easier to follow and more efficient. 

Subsequently, in step 512 of FIG. 8, the solution (u m+l , 
. . . , u) n and w^ is used in determining an optimal solution, 
w" to Problem Formulation |6.4| as generated according to 
I6.3J and Fact A.3. 

In step 516, the solution w n is used in defining the cost 
matrix C""""* 1 as in |6.7|. Subsequently, if n-m+l=2, then 
the assignment problem |6.8] may be solved straightfor- 
wardly using known techniques, such as forward/reverse 
auction algorithms. Following this, in step 528, the solution 
to the 2-dimensional assignment problem is assigned to the 
variable ASSIGNMT_SOLUTION and in step 532 
ASSIGNMT_SOLUTION is returned to a dimension three 
level recursion of MULTT_DIM_RELAX for solving a 
3-dimensional assignment problem. 

Alternatively, if in step 520, n-m+l>2, then in step 536 
the data structure(s) representing a problem formulation as 
in |6.8| is generated and assigned to the parameter, PROB_ 
FORMULATION. Subsequently, in step 540 a recursive 
copy of MULTLDIM _RELAX is invoked to solve the 
lower dimensional assignment problem represented by 
PROB_FORMULAXION. Upon the completion of step 
540, the parameter, ASSIGNMT__SOLUTION, contains the 
solution to the n-m+1 dimensional assignment problem. 
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Thus, in step 544. the n-m+1 solution is used to solve the 
n-dimensional assignment problem as discussed regarding 
equations |6.9|. Finally, in steps 548 and 552 the solution to 
the n-dimensional assignment problem is returned to the 

5 calling program so that, for example, it may be used in 
taking one or more actions such as (a) sending a warning to 
aircraft or sea facility; (b) controlling air traffic; (c) control- 
ling anti-aircraft or anti-missile equipment; (d) taking eva- 
sive action; or (e) surveilling an object. 

to n.3.1. Comments on the Various Procedures Provided by 
FIG. 8 

There are many procedures described by FIG. 8. One such 
procedure is the first embodiment of the multi-dimensional 
assignment solving.process of section 1. 1, That is, by speci- 

15 fying m=n-l in step 500, a single set of constraints is 
relaxed in step 508. TTius, one set of constraints incorporated 
is inta the objective function via the Lagrangian problem 
formulation, resulting in an m=n-i dimensional problem. 
The relaxed problem is subsequently maximized in step 512 

20 with respect to the corresponding Lagrange multipliers and 
then a feasible solution is reconstructed for the 
n-dimensional problem using a two-dimensional assignment 
problem. The second procedure provided by FIG. 8 is a 
novel approach which is not suggested by the first embodi- 

25 ment of section LI. In fact, the second procedure is some- 
what of a mirror image of the first embodiment in that n-2 
sets of constraints are simultaneously relaxed, yielding 
immediately an m=2 dimensional problem in step 512. Thus, 
a feasible solution to the n-dimensional problem is then 

30 recovered using a recursively obtained solution to an n-1 
dimensional problem via step 540. In this case, the function 
values and subgradients of ^^y 3 " optimally via a two 
dimensional assignment problem. The significant advantage 
here is. that there is no need for the merit or auxiliary 

35 functions, *¥ k embodiment of the multi-dimensional assign- 
ment solving process of section 1. 1 above and also there is 
no need for the more general merit or auxiliary functions 
H* subsection IL4.2 below. Further, note that all function 
values and subgradients used in the non-smooth maximiza- 

40 tion process are computed exactly (i.e., optimally) in this 
second procedure. Moreover, problem decomposition is now 
carried out for the n-dimensional problem; however, decom- 
position of the n- 1 dimensional recovery problem (and all 
lower recovery problems) is performed only after the prob- 

45 lem is formulated. 

Between these two procedures are a host of different 
relaxation schemes based on relaxing n-m sets of constraints 
to an m-dimensional problem (2<m<n), but these all have 
the same difficulties as the procedure for the first embodi- 

50 ment of section L 1 in that the relaxed problem is an NP-hard 
problem. To resolve this difficulty, we use an auxiliary or 
merit function ynotation m<n-l. the recovery 

procedure is still based on an NP-hard n-m+1 dimensional 
assignment problem. Note that the partitioning techniques 

55 similar to those discussed in Section 1.3 may be used to 
identify the assignment problem with a layered graph and 
then to find the disjoint components of this graph. In general, 
all relaxed problems can be decomposed prior to any non- 
smooth computations because their structure stays fixed 

60 throughout the-algorithm of FIG. 8. However, all recovery 
problems cannot be decomposed until they are formulated, 
as their structure changes as the solutions to the relaxed 
problems change. 

II.4. Details, and Refinements Relating to the Flowchart 

65 of FIG. 8 

Further explanation is provided here on how various steps 
of FIG. 8 are solved. Note that the refinements presented 
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here can significantly increase the speed of the relaxation 
procedure of step 508. 
D.4.1 Maximization of the Non-smooth Function 
One of the key steps in the Lagrangian relaxation algo- 
rithm in section IL3 is the solution of the problem 

Maximize {^V" 1 " *R**»1; - . . , n } [7.1] 

where u o *=0 for all k=m+l. . . . ji. The evaluation of 
(fr-y^ 1 "corresponding minimisa tion problem [6.2] Zn var- 
ies for each instance of (u""" 1 u"). The following 

discussion provides some properties of these functions. 

Fact A.4. Let u"* 1 " 1 u" be multiplier vectors associated 

with the (m+1)" through the n rt set of constraints [6.1]. let 
^mn Y 11 " objective function value of the n-dimensional 
assignment problem in equation |6.1]. let z n be any feasible 
solution of [6.1]. and let z" be an optimal solution of 16.1]. 



Then. O^y ' 



1 {u" 



. u"}and 



^Y" 1 "SV^SV.d") 



Furthermore. 



[7.2a] 



[7.2a] 



for all u*R Mt * with u o *=0 and k=m. . . . n. 

Most of the procedures for non-smooth optimization are 
based on generalized gradients called subgradients. given by 

the following definition. Definition. At us<u m+1 ? u") the 

set 64)^7 subdifferential of ^ 

lx. - . xR w ~l} [73] 

A vector qeS^^y 

If z" is an optimal solution of [6.2] computed during 
evaluation of O^y O^, £ yields the following i n -th com- 
ponent of a subgradient g of for O^y 



[7.4] 



*1 N k-l N k+l « 

for i k = 1, Nt and k =m+ 1, n. 



If z" is th e unique optimal solution of [6.2]. SO^y 
unique, then there are finitely many such solutions, say 
z"(l), . . . .z n (K). Given the corresponding subgradients, 

g l g* me subdifferential b<S>y {g l g*}. 

n.4.2 Mathematical Description of a Merit or Auxiliary 
Function 

For real-time needs, one must address the fact that the 
non-smooth optimization problem of step 508 (FIG. 8) 
requires the solution of an NP-hard problem for m>2. One 
approach to this problem is to use the following merit or 
auxiliary function to decide whether a function value has 
increased or decreased sufficiently in the line search or trust 
region methods: 



rrxS* 1 = 

dw*- 1 , ...V) 



[7.5J 



if m «= 2 

or (6.2) is soWcd optimally, 
, n") otherwise. 



where the multipliers u 3 u m that appear in lower order 

relaxations used to construct (suboptimal) solutions of the 
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m-dimensional relaxed problem (6.2) have been explicitly 
included. Note that y be solved optimally if m=2. For 
sufficiently small problems (6.2) or (6.4). one can more 
efficiently solve the NP-hard problem by branch and bound. 
5 This is the reason for the inclusion of the first case; 
otherwise, the relaxed function the merit function 
provides a lower bound for the optimal solution follows 
directly from Fact A.4 and the following fact 
Fact A.5. Given the definition of y 



17.6] 



The actual function value used in the optimization phase is 
4*mn y with the solution z (i n . . . i„ being a suboptimal solution 
constructed from a relaxation procedure applied to the 
15 m-dimensional problem. Again, the use of these lower order 
relaxed problems is the reason for the inclusion of the 
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multipliers u 

To explain how the merit function is used, suppose there 

is a current multiplier set (u 0U T* 1 u ot f ) and it is 

20 desirable to update to a new multiplier set (u„ 



via (u rt 



A " W.«" ouT 1 oJ* lot! .J"* obtained during the 
relaxation process used to compute a high quality solution to 
the relaxed m-dimensional assignment problem (6.2) at 

25 (u""" 1 u n )=(Uow mi " 1 u oid n ). The decision to accept 

(u^'ew. . . ji" Mw ) is then based on 4^ 3 old - old u^ 1 , 

• - • . Uoif^mnyne* ' neS ^AStODpittg Criteria 

commonly used in line searches. Again, (u^ u^J") 

represents the multiplier set used in the lower level relax- 
30 ation procedure to construct a high quality feasible solution 
to the m-dimensional relaxed problem (6.2) at (u^ 1 



Th e P oint is tnat toie one 



changes (u 1 ** 1 u") and uses the merit function W^y 3 

m /th-ijiA generally change the lower level multipliers (u 3 . 

35 U m ). 

An illustration of this merit function for m=n-l is given 
in "PARTITIONING MULTIPLE DATA SETS. MULTIDI- 
MENSIONAL ASSIGNMENTS AND LAGRANGIAN 
RELAXATION", by Aubrey B. Poore and Nenad Rijavec. 

40 and in "QUADRATIC ASSIGNMENT AND RELATED 
PROBLEMS. DIMACS SERIES IN DISCRETE MATH- 
EMATICS AND THEORETICAL COMPUTER 
SCIENCE", in American Mathematical Society. Providence. 
R.L. Vol. 16. 1994. pp. 25-37 edited by Panos M. Pardalos 

45 and Henry Wolkowicz. 

n.43 Non-smooth Optimization Methods 
By Fact A.4 the function ^y A affine. and concave, so 
that the negative of O^y A Thus the problem of maximizing 
O^y A optimization. Since there is a large literature on such 

50 problems, only a brief discussion of the primary classes of 
methods for solving such problems is provided. A first class 
of methods, known as subgradient methods, are reviewed 
and analyzed in the book by N. Z. Shor. Minimization 
Methods for Non-Differentiable Functions. Springer- Veriag, 

55 New York. 1985. Despite their relative simplicity, subgra- 
dient methods have some drawbacks that make them inap- 
propriate for the tracking problem They do not guarantee a 
descent direction at each iteration, they lack a clear stopping 
criterion, and exhibit poor convergence (less than linear). 

60 A more recent and sophisticated class of methods are the 
bundle methods; excellent developments are presented in the 
books of Hiriart-Urruty and Lemarechal: J.-B. Hiriart- 
Urruty and C. Lemarechal. Convex Analysis and Minimiza- 
tion Algorithms I and //. Springer- Veriag. Berlin. 1993; and 

65 the book of K. C. Kiwiel: Methods of Descent for Non- 
Differentiable Optimization, Lecture Notes in Mathematics* 
Vol. 1133. Springer-Veriag. Berlin. 1985. Bundle methods 
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retain a set (or bundle) of previously computed subgradients 
to determine the best possible descent direction at the 
current iteration. Because of the non-smoothness of 0^ 
may not provide a "sufficient" decrease in 0 WI bundle 
algorithms take a "null" step, wherein the bundle is enriched 
by a subgradient close to u k . As a result, bundle methods are 
non-ascent methods which satisfy the relaxed descent con- 
dition Q^y M A a+i*"** These methods have been 
shown to have good convergence properties. In particular, 
bundle method variants have been proven to converge in a 
finite number of steps for piece wise affine convex function- 
al in K.. Kiwiel, An Aggregate Subgradient Method for 
Non-smooth Convex Minimization, Mathematical Program- 
ming 27, (1983) pp. 320-341. and H. Schramm and J. Zowe, 
A version of the bundle idea for minimizing a non-smooth 
function: Conceptual idea, convergence analysis, numerical 
results. SIAM Journal on Optimization, 2 (1992), pp. 
121-152. 

All of the above non-smooth optimization methods 
require the computation A eS^^Y A the computation 
of the relaxed cost coefficients |6.3|. In both the first and 
second procedures discussed in section IL3. 1. maximization 
of 0^7 A In the most efficient implementations presently 
known of these two procedures, it was found that at least 
95% of the computational effort of the entire procedure is 
spent in the evaluation of the relaxed cost coefficients |6.3| 
as part of computing d> 2 „Y A that makes as efficient use of the 
subgradients and function values as is possible, even at the 
cost sophisticated manipulation of the subgradients. In 
evaluating three different bundle procedures: (a) the conju- 
gate subgradient method of Wolfe used in section 1.1 of the 
first embodiment of the present invention; (b) the aggregate 
subgradient method of Kiwiel (described in K. C. Kiwiel, An 
Aggregate Subgradient Method for Non-smooth Convex 
Minimization. Mathematical Programming 27, (1983) pp. 
320-341); and (c) the bundle trust method of Schramm and 
Zowe (described in H. Schramm and J. Zowe, A version of 
the bundle idea for minimizing a non-smooth function: 
Conceptual idea, convergence analysis, numerical results, 
SIAM Journal on Optimization. 2 (1992), pp. 121-152), it 
was determined that for a fixed number of non-smooth 
iterations, say, ten. the bundle-trust method provides good 
quality solutions with the fewest number of function and 
subgradient evaluations of all the methods, and is therefore 
the preferred method. 

n.4.4 The Two Dimensional Assignment Problem 

The forward/reverse auction algorithm of Bertsekas and 
Castanon (as described in D. P. Bertsekas and D. A. 
Castafion, A forward/reverse auction algorithm for asym- 
metric assignment problems. Computational Optimization 
and Applications, 1 (1992), pp. 277-298) is used to solve the 
many relaxed two dimensional problems that occur in the 
course of execution. 

n.4.5. Initial Multipliers and Hot Starts 

The effective use of "hot starts" is fundamental for 
real-time applications. A good initial set of multipliers can 
significantly reduce the number of non-smooth iterations 
(and hence the number of 0^ A quality recovered solution. 
Further, there are several ways that multipliers can be 
reused. First, if a n-dimensional problem is relaxed to a 
ra-dimensional problem, relaxation provides the multiplier 
set {u^V u"** 2 , . . . ,u"}. These can be used as the initial 
multipliers for the n-m+1 dimensional recovery problem for 
n-m+l>2. This approach has also worked well to reduce the 
number of non-smooth iterations during recovery. 

Further, for track maintenance, initial feasible solutions 
are generated as follows. When a new scan of information (a 
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new observation set) arrives from a sensor, one can obtain an 
initial primal feasible solution by matching new reports to 
existing tracks via a two dimensional assignment problem. 
This is known as the track -while -scan (TWS) approach. 

s Thus, an initial primal solution exists and then we use the 
above relaxation procedure is to make improvements to this 
TWS solution. Also for track maintenance, multipliers are 
available from a previously solved and closely related 
multi-dimensional assignment problems for all but the new 

10 observation set. 

n.4.6. Local Search Methods 

Given a feasible solution of the multi -dimensional assign- 
ment problem, one can consider local search procedures to 
improve this result, as described in C. H. Papadimitriou and 

15 K Stei glitz. Combinatorial Optimizption: Algorithms and 
Complexity. Prentice-Hall. Englewood Cliffs, N.J.. 1982, 
and J. Pearl, Heuristics: intelligent search strategies for 
computer problem solving, Addison- Wesley, 1984. One 
method is the idea of k-interchanges. Since for sparse 

20 problems only those active arcs that participate in the 
solution are stored, it is difficult to efficiently identify 
eligible variable exchange sets with some data structures for 
solving the assignment problem. However, a local search 
that may be more promising is to use the two dimensional 

25 assignment solver in the following way. Given a feasible 
solution to the multi-dimensional assignment problem, the 
indices that correspond to active arcs in the solution are 
enumerated. Subsequently, one coordinate is freed to for- 
mulate a two dimensional assignment problem with one 

30 index corresponding to the enumeration and the other to the 
freed coordinate, and then solve a two dimensional assign- 
ment problem to improve the freed index position. 
n.4.7. Problem Decomposition 

The procedures described thus far are all based on relax - 

35 ation. Due to the sparsity of assignment problems, however, 
frequently a decomposition of the problem into a collection 
of disjoint components can be done wherein each of the 
components can be solved independently. Due to the setup 
costs of Lagrangian relaxation, a branch and bound proce- 

40 dure is generally more efficient for small components, say 
ten to twenty feasible 0-1 variables (i.e., z i% J. Otherwise, 
the relaxation procedures described above are used. Perhaps 
the easiest way to view the decomposition method is to view 
the reports or measurements as a layered graph. A vertex is 

45 associated with each observation point, and an edge is 
allowed to connect two vertices only if the two observations 
belong to at least one feasible track of observations. Given 
this graph, the decomposition problem can then be posed as 
that of identifying the connected subcomponents of a graph 

50 which can be accomplished by constructing a spanning 
forest via a depth first search algorithm, as described in A. 
V. Aho, J. E. Hopcroft, and J. D. Ullman, The Design and 
Analysis of Computer Algorithms, Addison- Wesley Publish- 
ing Company, Reading, Mass., 1974. Additionally, decom- 

55 position of a different type might be based on the identifi- 
cation of bi- and tri-connected components of a graph and 
enumerating on the connections. Here is a technical expla- 
nation. Let 

Z={z iii2 i„tzili2 ... in is not pre assigned to zero} 

60 

denote the set of assignable variables. Define an undirected 
graph G(N*A) where the set of nodes is 

N=(z 4t "ln=l t ... 41; in=l Nm} 

65 and arcs, 

A^Cz^^jl^troil, j„*0, ji5<0 and there exists zili2 . . . ineZ such 
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, that and j,^}. 

Note that the nodes corresponding to zero index have not 
been Included in the above defined graph, since two vari- 
ables that have only the zero index in common can be 
assigned independently. Connected components of the graph 5 
are then easily found by constructing a spanning forest via 
a depth first search. (A detailed algorithm can be found in the 
book by Aho. Hopcroft and Ullman cited above). 
Furthermore, this procedure is used at each level in the 
relaxation, i.e.. is applied to each assignment problem |3.1] 10 
for n=3, . . . M. 

The original relaxation problem is decomposed first. All 
relaxed assignment problems can be decomposed a priori 
and all recover} 7 problems can be decomposed only after 
they are formulated. Hence, in the n-to-(n-l) case, we have is 
n-2 relaxed problems that can all be decomposed initially, 
and the recovery problems that are not decomposed (since 
they are all two dimensional). In the n-to-2 case, we have 
only one relaxed problem that can be decomposed initially. 
This case yields n-3 recovery problems, which can be 
decomposed only after they are formulated. 

The ever-increasing demand for sensor surveillance sys- 
tems is to accurately track and identify objects in real-time, 
even for dense target scenarios and in regions of high track 
contention. The use of multiple sensors, through more varied 
information, has the potential to greatly improve target state 
estimation and identification. However, to take full advan- 
tage of the available data, a multiple frame data association 
approach is needed. The most popular such approach is an 
enumerative technique called multiple hypothesis tracking 
(MHT). As an enumerative technique. MHT can be too 
computationally intensive for real-time needs because this 
problem is NP-hard. A promising approach is to utilize the 
recently developed Lagrangian relaxation algorithms for the 
multidimensional assignment problem; however, there are 
many other potentially good approaches to these assignment 
problems such as LP relaxation combined with an interior 
point method, GRASP, and parallelization. 

These data association problems are fundamentally 
important in tracking. Thus, our first objective is to present 
an overview of the tracking problem and the context within 
which these data association problems fit Following a brief 
review of the probabilistic framework for data association, 
the second objective is to formulate two models for track 
initiation and maintenance as multidimensional assignment 
problems. This formulation is based on a window moving 
over the frames of data. The first and simpler model uses the 
same window length for track initiation and maintenance, 
while the second model uses different lengths. As the 
window moves over the frames of data, one creates a 
sequence of these multidimensional assignment problems 
and there is an overlap region of frames common to adjacent 
windows. Thus, given the solution of one problem, one can 
develop a warm start for the solution to the next problem in 
the sequence, as shown hereinafter. Such information is 
critically important to the design of real-time algorithms. 

Overview of the Tracking Problem 

Tracking and data fusion are complex subjects of special- 
ization that can only be briefly summarized as they are 
related to the subject of this paper. The processing of track 
multiple targets using data from one or more sensors is 
typically partitioned into two stages or major functional 
blocks, namely, signal processing and data process. The first 
stage is the signal processing that takes raw sensor data and 
outputs reports. Reports are sometimes called observations, 
threshold exceedances. plots, hits, or returns, depending on 
the type of sensor. The true source of each report is usually 
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unknown and can be due to a target of interest, a false signal, 
or persistent background objects that can be moving in the 
field of view of the sensor. 

Two principal functions of data processing are tracking 
and discrimination or target identification. The discrimina- 
tion or target identification function distinguishes between 
tracks that are for targets of interest and other tracks such as 
those due to background. Also, if there is enough 
information, each track is classified as to the type of target 
it appears to represent. Target identification and cUscrimina- 
tion will not be discussed further in this paper except to 
comment here on the use of attributes (including identifica- 
tion information) in the data association. 

There are two types of attributes or features, namely, 
discrete valued variables and continuous valued variables. 
Available attributes of either or both types should be used in 
computing the log likelihoods or scores for data association. 
In discussing tracking in this paper, the attributes will not be 
mentioned explicitly but it should be understood that some 
reports and some tracks may include attributes information 
that should be and can be used in the track processing (both 
data association and filtering) if it is useful. 
Tracking Functions 

The tracking function accepts reports for each frame of 
data and constructs or maintains a target state estimate, the 
variance-covariance matrix of the state estimation error, and 
refined estimates (or probabilities) of the target attributes. 
The state estimate typically includes estimates of target 
position and velocity in three dimensions and possibly also 
30 acceleration and other variables of interest. 

The tracking function is typically viewed in two stages, 
namely, data association and filtering. Also, the process of 
constructing new tracks, called track initiation, is different 
from the process of updating existing tracks, called track 
35 maintenance. 

In track maintenance, the data association function 
decides how to relate the reports from the current frame of 
data to the previously computed tracks. In one approach, at 
most one report is assigned to each track, and in other 
40 approaches, weights are assigned to the pairings of reports 
to a track. After the data association, the filter function 
updates each target state estimate using the one or more 
(with weights) reports that were determined by the data 
association function. A filter commonly used for multiple 
45 target tracking and data fusion is the well know Kalman 
filter (or the extended Kalman filter in the case of nonlinear 
problems) or a simplified version of it 

In track initiation, typically a sequence of reports is 
selected (one from each of a few frames of data) to construct 
50 a new track. In track initiation, the filtering function con- 
structs a target state estimate and related information based 
on the selected sequence of reports. The new track is later 
updated by the track maintenance processing of the subse- 
quent frames of reports. 
55 In some trackers, there is also a tentative tracking function 
for processing recently established tracks until there is 
enough confidence to include them in track maintenance. 
While for simplicity of presentation, tentative tracking will 
not be included in the processing discussed in this paper, the 
60 techniques discussed could readily include one or more 
tentative tracking functions. 

There have been numerous approaches developed to 
perform the data association function. Since optimal data 
association is far too complex to implement, good but 
65 practical sub-optimal approaches are pursued. Data associa- 
tion approaches can be classified in a number of ways. On 
way to classify data association approaches is based on the 





5.959.574 



45 



46 



number of data frames used in the association process. In 
single frame data association for track maintenance, "hard 
decisions" are made each frame as to how the reports are to 
be related to the tracks. Some single frame approaches 
include: individual nearest neighbor. PDA, JPDA, and glo- 5 
bal nearest neighbor (sequential most probable hypothesis 
tracking), which uses a two dimensional assignment algo- 
rithm. In most single frame data association approaches, 
only one track per object is carried forward to be used in 
processing the next frame of reports. 10 

Multiple frame data association is more complex and 
frequently involves purposely carry forward more than one 
track per target to be used in processing the next frame of 
reports. By retaining more than one track per target, the 
tracking performance is improved at the cost of increased 15 
processing load. The best known multiple frame data asso- 
ciation approach is an enumerative technique called multiple 
hypothesis tracking (MHT) which enumerates all the global 
hypothesis with various pruning rules. 

As shown hereinafter, the log likelihood of the probability 20 
of a global hypothesis can be decomposed into the sum of 
scores of each track that is contained in the hypothesis. As 
a consequence, the most probable hypothesis can be iden- 
tified with the aid of a non-en umerative algorithm for the 
assignment so formulated. 25 

Interface Between Track Initiation And Track Mainte- 
nance 

There are a number of ways that track initiation and track 
maintenance functions can interact. Two ways that are 
especially pertinent to the methods of this paper will be 30 
discussed in this section. 

One approach is to treat these two functions sequentially, 
by first assigning reports to tracks and then using the 
information that remains to form new tracks. A better 
approach is to integrate both processes where initiating 35 
tracks compete for the new frame of reports on an equal 
basis, i.e., at the same time. 

In the integrated approach, the data association for track 
initiation and track maintenance processing are combined 
and conducted simultaneously. One assignment array is 40 
created that includes the track scores for potential tracks for 
both track initiation and track maintenance. The first dimen- 
sion of the assignment problem includes only all the tracks 
either created or updated in the processing of the prior 
frames of reports. Each of the remaining dimensions accom- 45 
modates one frame of reports. 

After the assignment algorithm finds a solution, each of 
the previously established tracks are updated with the report 
assigned to it from the second dimension of the assignment 
array. The remaining reports in the second dimension that 50 
were in the assignment solution are firmly established as 
track initiators. The unassigned reports in the second dimen- 
sion of the assignment array are discarded as false signals. 
The processing of the next frame of reports repeats this 
process using these updated tracks and the newly identified 55 
initiators in the first dimension of the cost array for process- 
ing the new frame of reports. Details of this approach are 
discussed hereinafter. 

The integrated approach just discussed, uses the same 
number of frames of reports for both track initiation and 60 
track maintenance. The goal in using the multi-dimensional 
assignment algorithm is to provide improved performance 
while minimizing the amount of processing required. 
. Typically, track initiation will benefit from more frames of 
reports than will track maintenance. Thus a second approach 65 
that integrates track initiation and track maintenance is 
discussed hereinafter, wherein the number of frames of 



reports is not the same for these two functions. This is a 
novel approach that is introduced in this paper. 
Multiple Sensor Processing 

There are many advantages and many ways to combining 
data from multiple sensors. There are also many ways of 
categorizing the different algorithmic architectures for pro- 
cessing data from multiple sensors. One approach outlines 
four generic algorithmic architectures. Two of these generic 
architectures are especially pertinent to this paper and are 
summarized briefly. 

In the Centralized Fusion algorithmic architecture, reports 
are combined from the various sensors to form global tracks. 
This algorithmic architecture is also called Central Level 
Tracking, Centralized Algorithmic Architecture, or simply 
the Type IV algorithmic architecture. In track maintenance 
for example, if a single frame data association is used, then 
for data association a frame of reports from one sensor is 
processed with the latest global tracks; then the global tracks 
are updated by the filter function. After the processing of this 
frame of reports is completed, a frame of the reports from 
another (or the same) sensor is processed with these updated 
global tracks. This process is continued as new frames of 
data become available to the system as a whole. 

In Centralized Fusion, using the multi-dimensional 
assignment algorithm for the data association with multiple 
sensors is similar to the processing of data from a single 
sensor. Instead of using multiple frames of reports from a 
single sensor, the multiple frames of reports come from 
multiple sensors. The frames of reports from all the sensors 
are ordered based on the nominal time each frame of reports 
is obtained and independent of which sensor provided the 
reports. In this way the approaches discussed in this paper 
can be extended to process reports from multiple sensors 
using the Centralized Fusion algorithmic architecture. 

The second pertinent algorithmic architecture is Track 
Fusion. This approach is also called the Hierarchical or 
Federated algorithmic architecture, sensor level tracking or 
simply the Type II algorithmic architecture. In track main- 
tenance for example, a processor for each sensor computes 
single sensor tracks; these tracks are then forwarded to the 
global tracker to compute global tracks based on data from 
all sensors. 

After the first time a sensor processor forwards tracks to 
the global tracker, then subsequent tracks for the same 
targets are cross-correlated with the existing global tracks. 
This track-to-track cross-correlation is due to the common 
history of the current sensor tracks and the tracks from the 
same sensor that were forwarded earlier to the global tracker. 
The processing must take this cross-correlation into account 
and there are a number of ways of compensating for this- 
cross-correlations. One method for dealing with this cross- 
correlation is to decorrelate the sensor tracks that are sent to 
the global tracker. There are a variety of ways to achieve this 
decorrelation and some are summarized in Drummond. 
1996, Signal Data Proc. of Small Targets 1995, 
-2561:369-383. 

Once the sensor tracks are decorrelated they can be 
processed by the global tracker in almost the same way as 
reports are processed. In the case of track fusion the asso- 
ciation process is referred to as track (or track-to-track) 
association rather than data (or report-to-track) association. 
If the sensor tracks are decorrelated, the global tracker can 
process the tracks from the various sensor processors in 
much the same way that the global tracker of Centralized 
Fusion processes reports. Accordingly the methods 
described in this paper can be readily extended to processing 
data from multiple sensors using either the Centralized 
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Fusion or Track Fusion or even a hybrid combination of both 
algorithmic architectures. 
Formulation of the Data Association Problem 
The goal of this section is to briefly outline the probabi- 
listic framework for the data association problems presented 
in this work. The technical details are presented elsewhere. 
The data association problems for multisensor and mill ti tar- 
get tracking considered in this work are generally posed as 
that of maximizing the posterior probability of the surveil- 
lance region (given the data) according to 
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-continued 

1 if ttij ) assigned to as a track, 

0 otherwise. 



Then, recognizing that 



10 



ftry|2* +l ) 



Maximize {/'(TVylZ"* 1 )^} 



(3.1) 



where Z**" 1 represents N+l data sets, y indices of the data 
(and thus induces a partition of the data). 



the problem (3.1) can be expressed as the N+l-dimensional 
assignment problem 



r\5 



15 



20 



Minimize 



-to+l w l"to+l 



Subject to y ... 2j -to+1 = M = 1 Wi. 



(3.5) 



partition, and P(TylZ* +1 ) is the posterior probability of a 
partition y N+l is defined below; however, this framework 
is currently sufficiently general to cover set packings and 25 
coverings. 

Consider N+l data sets Z(k)(k=l N+l). each 

consisting of M* actual reports and a dummy report Zq*. and 
let denote the cumulative data set defined by 

30 

2ik^{z ik k }»±0 and ZT+ l ={Z(l\ . . . , Z(*H)}, (3.2) 



to+l=° 



for it = 1, ... , M k and k = 2 N, 

^l- to +i 1" 1 -*- <*to+i 



€(0 t 1} for all it. 



respectively. (The dummy report serves several purposes 
in the representation of missing data, false reports, initiating 
tracks, and terminating tracks.) In multisensor data fusion 
and multitarget tracking the data sets Z(k) may represent 
different classes of objects, and each data set can arise from 
different sensors. For track initiation the objects are reports 
that must be partitioned into tracks and false alarms. In our 
formulation of track maintenance, which uses a moving 
window, one data set will be tracks and remaining data sets 
will be reports which are assigned to existing tracks, as false 
reports, or to initiating tracks. 

We specialize the problem to the case of set partitioning 
defined in the following way. Define a "track of data" as 

{z i( l ^i** 1 } wnere eacn h and z i* can assume the 

values of 0 and z©*. respectively. A partition of the data will 
refer to a collection of tracks of data wherein each report 
occurs exactly once in one of the tracks of data and such mat 
all data is used up; the occurrence of dummy report is 
unrestricted. The reference partition is that in which all 
reports are declared to be false. 

Next, under appropriate independence assumptions one 
can show 



where Cq , 0 is arbitrarily defined to be zero. Here, each 
group of sums in the constraints represents the fact that each 

35 non-dummy report occurs exactly once in a '"track of data.** 
One can modify this formulation to include multiassign- 
ments of one. some, or all the actual reports. The assignment 
problem (3.5) is changed accordingly. For example, if z* is 
to be assigned no more than, exactly, or no less than n, k 

40 times, then the "=!** in the constraint (3.5) is changed to "a. 



mLym ["J V-.to+i' 



<'i-to+i* 



a*- 1 is a likelihood ratio containing probabilities for 
detection, maneuvers, and termination as well as probability 
density functions for report errors, track initiation and ter- 
mination. Define 



«j »-to+i -to-r 



=, ^ n iA *," respectively. In making these changes, one must 
pay careful attention to the independence assumptions, 
which need not be valid in many applications. Expressions 
for the likelihood ratios can be found in the work of Poore, 
45 1994, Computational Optimization & AppL, 3:27-57. and 
the references therein. 
Track Initiation and Maintenance 
In this section we explain two multiframe assignment 
formulations to the track initiation and maintenance prob- 
50 lem. The continued use of all prior information is compu- 
tationally intensive for tracking, so that a window sliding 
over the frames of reports is used as the framework for track 
maintenance and track initiation within the window. The 
objectives are to describe an improved version of a simple 
55 method and then to put this into a more general framework 
(33) in which track initiation and maintenance have different 
length moving windows. 
The First Approach to Track Maintenance and Initiation 
The first method as explained in this section is an 
60 improved version of our first track maintenance scheme and 
uses the same window length for track initiation and main- 
tenance after the initialization step. The process is to start 
with a window of length N+l anchored at frame one. In the 
first step one there is only track initiation in that we assume 
(3 A) 65 no prior existing tracks. In the second and all subsequent 
frames, there is a window of length N anchored at frame k 
plus a collection of tracks up to frame k. This window is 
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denoted by {k; k+1 k+N}. The following explanation 

of the steps is much like mathematical induction in that we 
explain the first step and then step k to step k+1. 
Track Maintenance and Initiation: Step 1. Let 



(4.1) 



be an enumeration of all those zero-one variables in the 
solution of the assignment problem (3.5) (i.e., z i|i2 i^j) 
excluding all the false reports in the solution (i.e., all those 10 
zero-one variables with exactly one nonzero index) and 
zero-one variables in the solution for which (i 1 .i 2 )=(0-0). 
(The latter can correspond to tracks that initiate on frames 
three and higher.) These denote our initial tracks. 

Consider only the first two index sets in this enumeration 15 
(4.1) and add the zero index 1 2 =0 with the corresponding 
values of i A and i 2 being zero. Thus, the enumeration is now 
<ii(l 2 )a 2 (l 2 )}^0. The notation T^l^z.^.z^) will be 
used for this pairing. Suppose now that the next data set. i.e., 
the (N+2)"' set, is added to the problem. 2 o 

To explain the costs for the new problem, one starts with 
the hypothesis that a partition yeP* being true is now 
conditioned on the truth of the pairings on the first two 
frames being correct Corresponding to the sequence {T 2 
(l 2 ),z l3 z im2 ) , the likelihood function is then given by 25 



Ly= \~\ , t^ 3 ...* +2 where 



Next, define the cost and the corresponding zero-one vari- 
able 



CV3...^ + 2=- ln ^3-^ + 2' (4 ' 2b) 

1 if IzJ, z£,+\) is assigned to Track 7 2 </2), 

0 otherwise 



respectively. Then the track maintenance problem Maximize 
{LYryer* assignment 

Minimize > > ... ^ *V3-'w + 2* 2 '3 . '* + 2 



Subject to: > ... £ ^3.-^2 - 1. h = Uh ^ 

jpi 'w+2-° 

for i f »t M p and />«4 + 



50 

-continued 

Track Maintence and Initiation: Step k. At the beginning 
of the k" 1 step, we solve the following (N+l)-dimensional 
assignment problem. 



|4.4) 



Minimize: y y ... 



'i '1 + 1 • ■• k + w "t 't ♦ 1 - ■ • '* + N 



Subject to: \ . £ Vi + i -<i+* - 1. 4 - K - 

' K I. ... B D 



i 4+I -0 
*•! ^1+2 



/**0 ( 1+2 =0 'l+* =0 



^1+1 = 1, •■• Ml+I, 



(4.2 a) 



^fj-'w+a = Z ^3'"*w+2 " / ' 7 2"2>^3*"** + 2 * 30 



for i, = 1 M, and p = Jt+ W + Jk- 1, 

'1=0^+1-0 'i+i+yv-2- 0 



35 where for the first step 1, and L, are replaced by i t and 
respectively. The first index l k in the subscripts corresponds 
to the sequence of tracks {T 4 (y},/*=0 where T^Hz,/ 

(i t ) z^U} is a track of data from the solution of the 

problem prior to the formulation of (4.4) or if k=l then it is 

40 just the first data set in frame one. 

Suppose problem (4.4) has been solved and let the 
solution, i.e., those zero-one variables equal to one, be 
enumerated by 

««kiM*i('**i> W_ + i)))i w Ut - 1 < 4 ^ 

45 with the following exceptions. 

a. All zero one variables for which (l Jk i Jk+1 )=(0,0) are 
excluded. Thus, tracks that initiate on frames after the 
(k+1)" 1 are not included in the list. 

b. All zero-one variables whose subscripts have the form 
50 1*=0 and exactly one nonzero index in the remaining indices 

{ijt +1 , . . . ijt+yv} are excluded. These correspond to false 
reports on frames p=k+l k+N. 

c. All variables (Z/u. +1 *_) for which (i^,, . . . ♦ 

0) and 1 A)=6 m (4.5). In other words, the 

55 reports on the last N+l frames in the {T^l*),,^ 2 * 1 } are all 
dummy. Any solution with this feature corresponds to a 
terminated track. 

Given the enumeration (4.5). one now fixes the assignments 
on the first two index sets in the list (4.5). The zero index 
60 Ia+i^O IS added to the enumeration to specify (1^(0)4^1(0))= 
(0,0) that is used to represent false reports and tracks that 
initiate on frame k+2 or later, so that the enumeration (4.5) 
is now 

65 m^)JU^i))h M u ^ (46) 

Then, for 1^1=1, . . . X*+^ the \ tM th such track is denoted 
by {T t+1 (l A+l )={T ik (l Jk (l Jt+1 )),zk+l^ i (U 1 )}^and the (N+l>- 
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tuple iT^a^U^ 2 z, will denote a track 

T* +1 (W Plus a set of reports {z^ ^ h 

actual or dummy, that are feasible with the track T k +fik+i)- 



The (N+l)-tuple IWO).^ 



*} will 



denote a track that initiates in the sliding window. i.e.. on 5 
subsequent frames. A false report in the sliding window is 

one with all but one non-zero index ip for some p=k+2 

k+l+N in the (N+l)-tuple {T^(0).z Jt+2 * +2 . . . . 

The corresponding hypothesis about a partition y e T* 
being true is now conditioned on the truth of the L A+1 tracks 10 
existing at the beginning of the N-frame window. (Thus the 
assignments prior to this sliding window are fixed.) The 
likelihood function is given by 



Given this notation for the tracks and partition of the data in 
the frames {k-L . . . Jl . . . Je+J}. L T ^ ik) will denote the 
accumulated likelihood ratio up to and including frame 
p(p=k) for a track that is declared as existing on frame k as 
a solution of the assignment problem. In this notation, the 
likelihood for T ^(1*) and that of the association of 



v } with track T^y is given by 



'i+2 'i+l+N 



(4.7 a) 15 ^fSi^ViW ■ ■ ■ f0f my pSk_1, 

^K'i^i - - ■ W^WviW • ■ ■ WW ^ (4-10) 



20 



Next, define the cost and the zero-one variable by 

said at least one inert gas has an amount, on a molar basis, 
that is greater than the amount of said oxygen in said 
pressurized medium; 



(4.7b) 



1 if |< + + 2 2 i is assigned to T 1+ ,(/ t+1 y. 

0 otherwise 



respectively, so that the track extension problem, which was 
originally formulated as Maximize {LTflyer*}, can be 
expressed as exactly the same multi-dimensional assignment 
in (3.4) but with k replaced by k+1. Thus, we do not repeat 
it here. 

The Second Approach to Track Maintenance and Initia- 
tion 

In the first approach the window lengths were the same for 
both track maintenance and initiation. This can be inefficient 
in that one can usually use a shorter window for track 
maintenance than for track initiation. This section addresses 
such a formulation. 

The General Step k. To formulate a problem for track 
initiation and maintenance we consider a moving window 

centered as frame k of length I+J+ 1 denoted by [k-L 

k. . . . • k+J]. In this representation, the window length for 
track maintenance is J and that for initiation. I+J+l. The 
objective will be to explain the situation at this center and 
then the move to the same length window at center k+1 

denoted by [k+l-L . . . .k+1 Jc+l+J]. i.e.. by moving the 

frame one frame at time to the right. The explanation from 
the first step follows hereinafter. 

The notation for a track of data is 



(4J8) 



to 



respectively. 

The cost for the assignment of {z im * +1 . 
track Tj^y and the corresponding zero one variable are 
given by 



25 ^-W--^,^.^] (4 - U) 



30 



I: 



if -,$? N ] is assigned to track T k {l k \ 

otherwise 



respectively. Likewise, for costs associated with the false 
35 reports on the frames k-I to k and as associated with 
{z ihi * + \ . . . 3>i M ** J } a nd the corresponding zero-one vari- 
able are given by: 



40 



(4.12) 



1 * »**-/ *W> 

are assigned as a track, 
0 otherwise. 



45 



For the window of frames in the range [k-L ... i]. the 
easiest explanation of the partition of the reports is based on 
the definitions 



50 
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(4.13) 



i p I z? p belongs to one of the tracks listed on frame Jk, j 
it, to one of the T t (k) for some l k = 1, ... , Lt J 
/>jt*(U Af„l\r,jOU|0|, 



where the index \ k is used for an enumeration of those 
reports paired together. We also use the notation Tp^lJ to 
denote the sequence of reports belonging to track T^l*) but 
restricted to frames prior to and including p. Thus, 



so that all of those reports not used in the tracks T*(l*) on 
60 frame p are put into the set of available reports F pjc for the 
formation of tracks over the entire window. Thus, we 
formulate the assignment problem as 

Minimize 2 2 2 2^-/-'i4+ 1 -"W^-'-^ +I »^ 
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-continued decisions corresponding to the reports in our list of tracks 

«3 prior to and including frame k+l=I+2. This defines the k for 

subject to V* ... £ - ^+2 = ^ '2 = i. '2 problem (4.4) and one can then continue the development by 

rrl ^+2"° adding a frame to the window as in the general case. 



ij«0 



v^i "h+i If I+J>N, then one possibility is to start the process with 

/ k / /" ^'3 -»/ w+2 = 1. '3 = 1 M*. f,. ames an( j assuming J^N, proceed as before replac- 

<2=o ^-0 '" +2tt0 j n g 1 DV f or tne moment, and continue to add frames 

u *^ Wg-j *v» «„ +2 ' 1Q without lopping off the first frame in the window until 

/ / "/ / "' X *i'3-'*+2 " 1 reaches a window of length I+J+l. Then we proceed as in the 

/ 2 «o i 3 =o ip-i-oi p+1 «o previous paragraph. 

for «* = 1 A/ t and A = 2 /v, if I+J<N, then one can solve the track initiation problem 

th M 15 (3.5), formulate the problem with the center of the window 

V 1 ... V = *• " ■•* at k+1=N+1_ ^* enumerate the solutions as above, and lop off 

£0 ^+1-1-0 the first N-J-I frames. Then, we proceed just as in the case 



I+J=N. 



l/v+i € |0, 1| for all i'i i/v+i 



20 A primary objective in this work has been to demonstrate 
how multidimensional assignment problems arise in the 

Suppose problem (4-14) has been solved and let the tR| environment. The problem of track initiation and 

solution. i.e.. .hose zero-one vanables equal to one. be maintenance has been formulated within the framework of a 

l (4 15) moving window over the frames of data. The solution of 

/ T«™«*ed by \ *+i these NP-hard, noisy, large scale, and sparse problems to the 

4+1 noise level in the problem is fundamental to superior track 

Izf . , . . u , . . A^+i estimation and identification. Thus, one must utilize the 

special structure in the problems as well take advantage of 
special information that is available. Since these moving 

with the following exceptions. 30 windows are overlapping, there are some algorithm efficien- 

a. All zero one variables in the second list for which cies that can be identified and that take advantages of the 

(h-i* • • • J A+iMO 0,0) are excluded. Thus, tracks that overlap in the windows from one frame of reports to the 

initiate on frames after the (k+1)'* are not included in the list. nex t. Here is an example of the use of a primal solution of 

b. All false reports axe excluded, i.e., all zero-one van- one pro blem to warm start the solution of the next problem 
ables in the second list whose subscripts have exactly one 35 - n ^ se q Uence 

nonzero index. 

c. All variables z,* i T for which (i^,, . . . a* + *)= Suppose we have solved problem (4.4) and have enumer- 

(0,0 0) and z^nT^l*)^ for p=k-K, . . . ,k where K^0 ated ^ those zero-one variables in the solution of (4.4) as 

is user specified. Thus the track T^IJ is terminated if it is in (4>5 ) Add tne ^ m dex 1*^=0, so that the enumeration 

not observed over K+J+l frames. 40 | s 

Given the enumeration (4.15). one now fixes the assign- 
ments on the all index sets up to and including the (k+l) th 

index sets. {(W^iW.Oui) ^mOm))} (5-0 

( T *l'*<'^>.4/ + V'*+»>) (416> 45 With this enumeration one can define the cost by 

for = 1 li, 

«;> + .x .... <<4*ix w^<ww . . • JMnvM* < 52 > 



for U+i = 1 *4+i + 1 i-i+i 



50 



and the two dimensional assignment problem 



Thus one can now formulate the assignment problem for 

the next problem exactly as in (4-14) but with k replaced by < 5 - 3 > 

k+i. Thus, we do not repeat it here. ^ 3 ^ 2^ h A+i'i+uw c Vj <^ 

The Initial Step. Here is one explanation for the initial 55 'a+i-°a+i+* 

step. First, assume that N=I+J. In this case, we start the track ^ M ^ + * 
initiation with a solution of (3.5). Let 



Subject to £ ^^-W^-l, 



be an enumeration of the solution set of (3.5), i.e., those 60 £ 'Wi+i** a Xt ° 1 

zero-one variables z^m^) , . . wd) =1 ' including Zqq <f=l '*+«"° 

corresponding to 1 2 =0, but excluding all those zero-one ^ ( e|0, i| forali // 

variables that are assigned to one and correspond to false **i 

reports (i.e., there is exactly one nonzero index in the 

subscript of z i(i2 _ t 4 ), all those zero-one variables that are 65 

assigned to one andTcorrespond to tracks that initiate on Let w be an optimal or feasible solution to this two- 
frames higher than 1+2. Then we fix the data association dimensional assignment problem and define 
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1 if 



1 for some = 1, L*+i 



or if (Wt + i+w) = O,0 
0 otherwise. 



This need not satisfy the constraints in that there are usually 
many objects left unassigned. Thus, one can complete the 
assignment by using the zero-one variables in (4.4) with k 
replaced by k+1 with exactly one nonzero index correspond- 
ing to any unassigned object or data report 

For the dual solutions, the multipliers arising from the 
solution of the two dimensional assignment problem (5.3) 
corresponding to the second variable, i.e.. { u i fcfW / M " 1+A '}< 4tH 
n Mm +"=G. These are good initial values for use in a relaxation 
scheme [11.12]. Finally, note that one can also develop a 
warm start for problem (4.14) in a similar fashion. 

Multidimensional assignment problems govern the cen- 
tral problem of data association in multisensor and multi- 
target tracking. i.e.. the problem of partitioning observations 
from multiple scans of the surveillance region and from 
either single or multiple sensors into tracks and false alarms. 
This fundamental problem can be stated as 



Minimize 



Subject to: \ ... £ V*'* = ^i-M-™*! • 

Z-Z Z - ^*"^^ 1 

for i p = 1, ...,M, and /><= 2, 1, 

li v -.i N e{0, 1) for all ii i*. 

for i t = 1 Mt and k = 2, ... , N, 

2i r ..liv+i € 10, 1| for all ii 



methods for solving this NP-hard problem optimally are 
(5.4) enumerative in nature, with branch -and-bound being the 
most efficient; however, such methods are much too slow for 
real-time applications. Thus one must resort to suboptimal 
5 approaches. Ideally, any such algorithm should solve the 
problem to within the noise level, assuming, of course, that 
one can measure this noise level in the physical problem and 
the mathematical method provides a way to decide if the 
criterion has been met 
10 There are many algorithms that can be used to construct 
suboptimal solutions to NP-hard combinatorial optimization 
problems. These include greedy (and its many variants), 
relaxation, simulated annealing, tabu search, genetic 
algorithms, and neural netork algorithms. For the three 
15 dimensional assignment problem. Pierskalla developed the 
tri-substitution method, which is a variant of the simplex 
method. Frieze and Yadegar introduced a method based on 
Lagrangian relaxation in which a feasible solution is recov- 
ered using information provided by the relaxed solution. Our 
20 choice of approaches is strongly influenced by the need to 
balance real-time performance and solution quality. 
Lagrangian relaxation based methods have been used suc- 
cessfully in prior tracking applications. An advantage of 
these methods is that they provide both an upper and lower 
25 bound on the optimal solution, which can then be used to 
measure the solution quality. These works extend the 
method of Frieze and Yadegar to the multidimensional case. 
Probabilistic Framework for Data Association. (ABP) 
(i.D The goal of this section is to explain the formulation of the 
30 data association problems that governs large classes of data 
association problems in centralized or hybrid centralized- 
sensor level multisensor/multitarget tracking. The presenta- 
tion is brief; technical details are presented for both track 
initiation and maintenance in the work of Poore. The for- 
35 mulation is of sufficient generality to cover the MHT work 
of Reid. Blackman and Stein and modifications by Kurien to 
include maneuvering targets. As suggested by Blackman. 
this formulation can also be modified to include target 
features (e.g.. size and type) into the scoring function. The 
40 recent work of Poore and Drummond significantly extends 
the work of this section to new approaches for multiple 
sensor centralized tracking. Future work will involve exten- 
sions to track-to-track correlation. 
The data association problems for multisensor and mul- 
45 titarget tracking considered in this work are generally posed 
as that of maximizing the posterior probability of the sur- 
veillance region (given the data) according to 



Maximize {Pfr^iif W} 



(2.1) 



where Cq 0 is arbitrarily defined to be zero and is included 50 
for notational convenience. One can modify this formulation 
to include multiassignments of one. some, or all of the actual 
reports. The zero index is used in representing missing data, 
false alarms, initiating and terminating tracks. In these 
problems, we assume that all zero-one variables z if _i N with 55 
precisely one nonzero index are free to be assigned and that 
the corresponding cost coefficients are well-defined. (This is 
a valid assumption in the tracking environment) Although 
not required, these cost coefficients with exactly one nonzero 
index can be translated to zero by cost shifting without 60 
changing the optimal assignment. Finally, this formulation is 
of sufficient generality to include the symmetric problem 
and the asymmetric inequality problem. 

The data association problems in tracking that are formu- 
lated as (1.1) have several characteristics. They are normally 65 
sparse, the cost coefficients are noisy and the problem is 
NP-hard, but it must be solved in real-time. The only known 



where Z N represents N data sets, y of the data (and thus 
induces a partition of the data). F6 



F(T ylZ") is the posterior probability of a partition y true 
given the data Z*. The term partition is defined below; 
however, this framework is currently sufficiently general to 
cover set packings and coverings. 

Consider N data sets Z(k) (k=L . . . .N) each of M* reports 
{Zi k } tk Mk =l* and let Z N denote the cumulative data set 
deAned by 



(2a) 
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^ s n l » ~ >«' 



■4 



Minimize £ c 'i-'**i-i* 

Subject to ^ = 1 *'» = 1. ••• » M|) 

2) -1 ('2 - 1 M 2 ) 

M, and p = 2,.../V-l) 

Z D 1 ° 1 M *>" 

e{0, 1| for aU /, i N . 
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respectively. In multisensor data fusion and multitarget 
tracking the data sets Z(k) may represent different classes of 
objects, and each data set can arise from different sensors. 
For track initiation the objects are measurements that must 
be partitioned into tracks and false alarms. In our formula- 
tion of track maintenance, which uses a moving window 
over time, one data set will be tracks and remaining data sets 
will be measurements which are assigned to existing tracks, 
as false measurements, or to initiating tracks. In sensor level 
tracking, the objects to be fused are tracks. In centralized 
fusion, the objects may all be measurements that represent 
targets or false reports, and the problem is to determine 
which measurements emanate from a common source. 

We specialize the problem to the case of set partitioning 
defined in the following way. First, for notational conve- 
nience in representing tracks, we add a zero index to each of 
the index sets in (2.2). a dummy report to each of the data 
sets Z(k) in (2.2). and define a "track of data" as (z./.z,^) 
where i k and z£ can now assume the values of 0 and z^ 
respectively. A partition of the data will refer to a collection 20 
of tracks of data wherein each report occurs exactly once in 
one of the tracks of data and such that all data is used up; the 
occurrence of dummy report is unrestricted The dummy 
report Zq* serves several purposes in the representation of 
missing data, false reports, initiating tracks, and terminating 25 
tracks. The reference partition is that in which all reports are 
declared to be false. 

Next under appropriate independence assumptions one 
can show 



10 



15 
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(2.3) 



L (i ijv is a likelihood ratio containing probabilities for 35 
detection, maneuvers, and termination as well as probability 
density functions for measurement errors, track initiation 
and termination. Then if c^ iw =-lnL J( . . . ^ 



signments of one. some, or all the actual reports. The 
assignment problem (2.5) is changed accordingly. For 
example, if z ik k is to be assigned no more than, exactly, or no 
less than n, 4 A times, then the "=1" in the constraint (2.5) is 
changed to ^, =, ^n <A *" respectively. Modifications for 
group tracking and multiresolution features of the surveil- 
lance region will be addressed in future work. In making 
these changes, one must pay careful attention to the inde- 
pendence assumptions that need not be valid in many 
applications. 

Expressions for the likelihood ratios v can be 

found in the work of Poore. These expressions include the 
developments of Reid, Kurien, and Stein and Blackman. 
What's more, they are easily modified to include target 
features and to account for different sensor types. In track 
initiation, the N data sets all represent reports from N 
sensors, possibly all the same. For track maintenance, we 
use a sliding window of N data sets and one data set 
containing established tracks. The formulation is the same as 
above except that the dimension of the assignment problem 
is now N+l.@@@ 

To explain the costs for the new problem, one starts with 
the hypothesis that a partition y e T* conditioned on the truth 
of the pairings on the first two frames being correct. Cor- 
responding to the sequence {T 2 (l 2 ), z iy z,^ 2 }< the 

likelihood function is then given by 

r>= j~j L, 2h .- UlM where 

Lt 2 =1. 



Next, define the cost and the corresponding zero-one vari- 
able 



Using (2.3) and the zero-one variable z f| t =l if (i,, . . . , 
i^)e y and 0 otherwise, one can then write the problem (2. 1) 45 
as the following N-dimensional assignment problem: 



(2.4) 40 C^-w^^-^...^. 

1 if I* <&» * 

z *2h-itt*i = assigned to Track T 2 (/ 2 ), 

0 otherwise 



(4.2 b) 



(2.5) 



50 



respectively. Then the track maintenance problem Maximize 
{Ijyryer* assignment 

/ 2 -0 
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Subject to: } — 2j ***** "■'*+> ° ' 2 ° l * ' 2 



where Cq 0 is arbitrarily defined to be zero. Here, each 
group of sums in the constraints represents the fact that 65 
each-non-dummy report occurs exactly once in a 'track of 
data.** One can modify this formulation to include multias- 



for i k - 1 M k and 2,... N t 
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-continued 



60 



Then, for L + .=l 



th 



w (^iHTM w ))V,l. 



by T, 

{T^iv^i7**-, w .... 
T w (l*».i) plus a set of reports {z, 



such track is denoted 
and the (N+l)-tuple 
will denote a track 



5 actual or dummy, that are feasible with the track T^O^). 



2fj...lyv+i € {0, II for all i'i 



The (N+l)-tuple <T W <0). z,- 



v }. will 



Track Maintence and Initiation: Step k. At the beginning 
of the k** step, we solve the following (N+l>dimensional 
assignment problem- 



denote a track that initiates in the sliding window, i.e.. on 
subsequent frames. A false report in the sliding window is 
one with all but one non-zero index L for some p=k+2. 



10 k+l+N in the (N+l)-tuple {T^,(0). z, 



Minimize 



(4.4) 



Subject to: 



' 2 



U 4 - 1 U, 



The corresponding hypothesis about a partition y e T* 
being true is now conditioned on the truth of the L* +1 tracks 
existing at the beginning of the N-frame window. (Thus the 
assignments prior to this sliding window are fixed.) The 
15 likelihood function is given by 



0 *k+N m 



where 



(4.7 a) 



Zj Zj " sl - 

= 1, Mt+L 
'1-0 '* + * 

for i, = 1 M p and /> = k+ 2, ... , N, +k + 1, 

i 4 «0/ t+I «0 '4+1+^-2"° 

° 1 

• W € H for all / A , i'a+jv 
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J i+2 U+l+W 

tr t+ , ( 0) = 1. 

Next, define the cost and the zero-one variable by 



1 if l£* 



't+2* 

0 otherwise. 



(4.7 b) 
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where for the first step l x and L A are replaced by ij and M 2 . 
respectively. The first index 1* in the subscripts corresponds 
to the sequence of tracks {T A (ljt) },/*=() where T^lJ^z,/ 

(ljt) 2,^(1*)} is a track of data from the solution of the 

problem prior to the formulation of (4.4) or if k=l then it is 
just the first data set in frame one. 

Suppose problem (4.4) has been solved and let the 
solution. i.e.. those zero-one variables equal to one, be 
enumerated by 

with the following exceptions. 

a. All zero one variables for which (l^i^MO-O) are 
excluded. Thus, tracks that initiate on frames after the 
(k+l) r/l are not included in the list. 

b. All zero-one variables whose subscripts have the form 
1 A =0 and exactly one nonzero index in the remaining indices 
{i*^, . . . l k + N } are excluded. These correspond to false 
reports on frames p=k+L . . . Jc+N. 



respectively, so that the track extension problem, which was 
originally formulated as Maximize {L^ly^S as exactly the 
same multi-dimensional assignment in (3.4) but with k 
replaced by k+1. Thus, we do not repeat it here. 

The Second Approach to Track Maintenance and Initia- 
tion 

In the first approach the window lengths were the same for 
bom track maintenance and initiation. This can be inefficient 
in that one can usually use a shorter window for track 
maintenance than for track initiation. This section addresses 
such a formulation. 

The General Step L To formulate a problem for track 
initiation and maintenance we consider a moving window 
centered as frame k of length I+J+l denoted by [k-L .... 
k, . . . Jc+J]. In mis representation, the window length for 
track maintenance is J and that for initiation. I+J+l. The 
objective will be to explain the situation at this center and 
then the move to the same length window at center k+1 

denoted by [k+l-L . . . Jc+1 Jc+l+J]. i.e.. by moving the 

50 frame one frame at time to the right. The explanation from 
the first step follows hereinafter. 
The notation for a track of data is 
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(4.8) 



c. All variables z, 



a^N) for which (ik+1. 



55 



ik+NM0,0 0) and i^l^H* in (4.5). In other words, the 

reports on the last N+l frames in the {T^IJ. z^*"* 1 } are all 
dummy. Any solution with this feature corresponds to a 
terminated track. 

Given the enumeration (4.5). one now fixes the assignments 
on the first two index sets in the list (4.5). The zero index 
1^=0 is added to the enumeration to specify M0)1 M (0))= 
(0.0) that is used to represent false reports and tracks that 
initiate on frame k+2 or later, so that the enumeration (4.5) 
is now 

i(tt/*iV*i(ki)>Wt^« (4.6) 



where the index 1* is used for an enumeration of those 
reports paired together. We also use the notation T^l*) to 
denote the sequence of reports belonging to track T*(l J but 
restricted to frames prior to and including p. Thus. 

n&H vow* • • ■ V<'*» 

nj£l k >T P 4l k yj{z^o<r l Wtf*" 1 W> fOT pSk-l.(«) 

65 Given this notation for the tracks and partition of the data in 

the frames {k-I Jc . . . Jc+J}. ^TpjtOjt) will denote the 

accumulated likelihood ratio up to and including frame 
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p(p^k) for a track that is declared as existing on frame k as 
a solution of the assignment problem. In this notation, the 
likelihood for T,,^) and that of the association of 



(4.10) 



. . & iMI \ with track T k (\ k ) is given by 
L nwr L W k ^ li i L \ . . . naa for any p£k-l, 

respectively. 

The cost for the assignment of {z iitl * +1 2 iw***} to 

track T k (l k ) and the corresponding zero one variable are 
given by 



10 



1 1< + -', <SI 



0 otherwise. 



14.11) 15 
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-continued 



Subject ^ Z ZZC-V, 



for '-jei^jtMOl and $ = ...jfc. 



' 1, ... , m a 



and $ = jfc + 1, ...,* + 7, 



'i. U+i. 



for $ = A + 1 Jt + y. 



is assigned to track T^ttJ, respectively. Likewise, for costs 
associated with the false reports on the frames k-I to k and 



as associated with { z, 



'iti 



7 * 



ing zero-one variable are given by: 

1 l if V/ fti-;-- 



} and the correspond- 



Suppose problem (4.14) has been solved and let the 
25 solution, i.e., those zero-one variables equal to one, be 
enumerated by 



(4.12) 



(4.15) 



30 K-i^+^-'it'i+iWi+ifi+^-Wfi+iili** 1 r 



are assigned as a track, 
0 otherwise. 



For the window of frames in the range |k~I, . . . ,k], the 
easiest explanation of the partition of the reports is based on 
the definitions 



(4.13) 



belongs to one of the tracks listed on frame it, 
i.e., to one of the T k (l k ) for some l k = 1 1+ 

f,jt=(|l M P \\T pJt )U[0l 



with the following exceptions. 

35 

a. All zero one variables in the second list for which 
(i w ^U+iMO* • • ■ .0,0) are excluded. Thus, tracks that 
initiate on frames after the (k+1)" 1 are not included in the list. 

40 b. All false reports are excluded, i.e., all zero-one vari- 
ables in the second list whose subscripts have exactly one 
nonzero index. 



45 



so that all of those reports not used in the tracks 1^(1*) on 
frame p are put into the set of available reports for the 
formation of tracks over the entire window. Thus, we 
formulate the assignment problem as 



c. All variables z Vw . . . , for which (i k+1 W<0. 

0 0) and Z(p)nT A (U={0} for p=k-K, . . . ,k where K£0 

is user specified. Thus the track T^l*) is terminated if it is 
not observed over K+J+l frames. 

Given the enumeration (4.15), one now fixes the assign- 
ments on the all index sets up to and including the (k+1)'* 
index sets. 



(4.16) 



Mil 



Thus one can now formulate the assignment problem for 
the next problem exactly as in (4. 14) but with k replaced by 
60 k+1. Thus, we do not repeat it here. 
Zj i 2j 2j Zi^-/- ■'W* + r-W z <fc-/ "Vi+r-'*w + The Initial Step. Here is one explanation for the initial 

p pa n,i+1 stC p t assume that N=I+J. In this case, we start the track 

initiation with a solution of (3.5). Let 

Z Z S^i-w^i-w 65 WkWh) • • • MW^Mh*** 

i 4 - In .i+iir»o ^ an enumeration of the solution set of (3.5), Le., those 

zero-one variables z iiih)hil0 . . . Wfl) =l, including z^ 0 =1 
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corresponding to 1 2=0. but excluding all those zero-one 
variables that are assigned to one and correspond to false 
reports (i.e.. there is exactly one nonzero index in the 
subscript of z iih i„+ iy all those zero-one variables that are 
assigned to one and correspond to tracks that initiate on 
frames higher than 1+2. Then we fix the data association 
decisions corresponding to the reports in our list of tracks 
prior to and including frame k+ 1=1+2. This defines the k for 
problem (4.4) and one can then continue the development by 
adding a frame to the window as in the general case. 

If I+J>N. then one possibility is to start the process with 
N+l frames, and assuming J=N. proceed as before replac- 
ing I by N-J for the moment, and continue to add frames 
without lopping off the first frame in the window until 
reaches a window of length I+J+l. Then we proceed as in the 
previous paragraph. 

If I+J<N. then one can solve the track initiation problem 
(3.5). formulate the problem with the center of the window 
at k+l=N+l-J. enumerate the solutions as above, and lop off 
the first N-J-I frames. Then, we proceed just as in the case 
I+J=N. 

A primary objective in this work has been to demonstrate 
how multidimensional assignment problems arise in the 
tracking environment. The problem of track initiation and 
maintenance has been formulated within the framework of a 
moving window over the frames of data. The solution of 
these NP-hard. noisy, large scale, and sparse problems to the 
noise level in the problem is fundamental to superior track 
estimation and identification. Thus, one must utilize the 
special structure in the problems as well take advantage of 
special information that is available. Since these moving 
windows are overlapping, there are some algorithm efficien- 
cies that can been identified and that take advantages of the 
overlap in the windows from one frame of reports to the 
next Here is an example of the use of a primal solution of 
one problem to warm start the solution of the next problem 
in the sequence. 

Suppose we have solved problem (4.4) and have enumer- 
ated all those zero-one variables in the solution of (4.4) as 
in (4.5). Add the zero index ljt +1 =0. so that the enumeration 
is 



With this enumeration one can define the cost by 

and the two dimensional assignment problem 
fl^ B Minimis £ Z ^iW^iW 3 ^ 



(5-1) 



(5-2) 
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2l 



(5.4) 



1 if W = '*+*('*+! H 

and H>/ 1+J / t+i+A , - for some l M = 1, 1 

or if (/*+!, ii + i+w) = <0. Oi, 

0 otherwise. 
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This need not satisfy the constraints in that there are 
usually many objects left unas signed. Thus, one can com- 
plete the assignment by using the zero-one variables in (4.4) 
15 with k replaced by k+1 with exactly one nonzero index 
corresponding to any unassigned object or data report 

For the dual solutions, the multipliers arising from the 
solution of the two dimensional assignment problem (5.3) 
corresponding to the second variable, i.e.. \ u i M J^ L ^ii M ^ 
20 ^Af^HN^Q These are good initial values for use in a relaxation 
scheme [11.12]. Finally, note that one can also develop a 
warm start for problem (4.14) in a similar fashion. 

Multidimensional assignment problems govern the cen- 
25 tral problem of data association in multisensor and multi- 
target tracking, i.e.. the problem of partitioning observations 
from multiple scans of the surveillance region and from 
either single or multiple sensors into tracks and false alarms. 
This fundamental problem can be stated as 
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Minimize 



f 1 v 

Subject to: \ £ l^...^ 



(1.1) 



1. i,-!.!...^! 



*2«0 

for i p = 1, ... , M p and p- 2, ... » N - 1, 

y z *r = i - iN = 1 M *' 

Ziy»i N €(0,1} for all i t i*. 

for i t = 1 M k and *=2,...,7V, 

* r ..l A+1 €{0,1| for all i, 
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Let w be an optimal or feasible solution to this 
dimensional assignment problem and define 



where c 0 0 is arbitrarily defined to be zero and is included 
(5 J) for notational convenience. One can modify this formulation 
55 to include multiassignments of one. some, or all of the actual 
reports. The zero index is used in representing missing data, 
false alarms, initiating and terminating tracks. In these 
problems, we assume that all zero-one variables z j; i^ with 
precisely one nonzero index are free to be assigned and that 
the corresponding cost coefficients are well-defined. (This is 
a valid assumption in the tracking environment.) Although 
not required, these cost coefficients with exactly one nonzero 
index can be translated to zero by cost shifting without 
65 changing the optimal assignment Finally, this formulation is 
two- of sufficient generality to include the symmetric problem 
and the asymmetric inequality problem. 
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The data association problems in tracking that are formu- 
lated as (1.1) have several characteristics. They are normally 
sparse, the cost coefficients are noisy and the problem is 
NP-hard, but it must be solved in real-time. The only known 
methods for solving this NP-hard problem optimally are 
enumerative in nature, with branch-and-bound being the 
most efficient; however, such methods are much too slow for 
real-time applications. Thus one must resort to suboptimal 
approaches. Ideally, any such algorithm should solve the 
problem to within the noise level, assuming, of course, that 
one can measure this noise level in the physical problem and 
the mathematical method provides a way to decide if the 
criterion has been met. 

There are many algorithms that can be used to construct 
suboptimal solutions to NP-hard combinatorial optimization 
problems. These include greedy (and its many variants), 
relaxation, simulated annealing, tabu search, genetic 
algorithms, and neural netork algorithms. For the three 
dimensional assignment problem. Pierskalla developed the 
tri-substitution method, which is a variant of the simplex 
method. Frieze and Yadegar introduced a method based on 
Lagrangian relaxation in which a feasible solution is recov- 
ered using information provided by the relaxed solution. Our 
choice of approaches is strongly influenced by the need to 
balance real-time performance and solution quality. 
Lagrangian relaxation based methods have been used suc- 
cessfully in prior tracking applications. An advantage of 
these methods is that they provide both an upper and lower 
bound on the optimal solution, which can then be used to 
measure the solution quality. These works extend the 
method of Frieze and Yadegar to the multidimensional case. 

Probabilistic Framework for Data Association. (ABP) 

The goal of this section is to explain the formulation of the 
data association problems that governs large classes of data 
association problems in centralized or hybrid centralized- 
sensor level multisensor/multitarget tracking. The presenta- 
tion is brief; technical details are presented for both track 
initiation and maintenance in the work of Poore. The for- 
mulation is of sufficient generality to cover the MHT work 
of Reid, Blackman and Stein and modifications by Kurien to 
include maneuvering targets. As suggested by Blackman, 
this formulation can also be modified to include target 
features (e.g., size and type) into the scoring function. The 
recent work of Poore and Drummond significantly extends 
the work of this section to new approaches for multiple 
sensor centralized tracking. Future work will involve exten- 
sions to track-to-track correlation. 

The data association problems for multisensor and mul- 
titarget tracking considered in this work are generally posed 
as that of maximizing the posterior probability of the sur- 
veillance region (given the data) according to 



Maximize </>(T=ylZ*tyy«r* } 
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Consider N data sets Z(k) (k=I .N) each of M k reports 

{z.^.^1, and let Z" denote the cumulative data set 
defined by 



ZW=Ui/hA*" =I **** ^ZO), - * • 



(2-2) 



respectively. In multisensor data fusion and multitarget 
tracking the data sets Z(k) may represent different classes of 
objects, and each data set can arise from different sensors. 
For track initiation the objects are measurements that must 
be partitioned into tracks and false alarms. In our formula- 
tion of track maintenance, which uses a moving window 
over time, one data set will be tracks and remaining data sets 
will be measurements which are assigned to existing tracks, 
as false measurements, or to initiating tracks. In sensor level 
tracking, the objects to be fused are tracks. In centralized 
fusion, the objects may all be measurements that represent 
targets or false reports, and the problem is to determine 
which measurements emanate from a common source. 

We specialize the problem to the case of set partitioning 
defined in the following way. First, for notational conve- 
nience in representing tracks, we add a zero index to each of 
the index sets in (2.2), a dummy report to each of the data 
sets Z(k) in (2.2), and define a 'track of data" as (z, Vz,/) 
where @@i A and z,.* can now assume the values of 0 and 
z 0 *. respectively. A partition of the data will refer to a 
collection of tracks of data wherein each report occurs 
exactly once in one of the tracks of data and such that all data 
is used up; the occurrence of dummy report is unrestricted. 
The dummy report Zq* serves several purposes in the rep- 
resentation of missing data, false reports, initiating tracks, 
and terminating tracks. The reference partition is that in 
which all reports are declared to be false. 

Next under appropriate independence assumptions one 
can show 



P[V= y°\Z») 



(2.3) 
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L (| iyv is a likelihood ratio containing probabilities for 
detection, maneuvers, and termination as well as probability 
density functions for measurement errors, track initiation 



45 and termination. Then if c, i^^JnL, 



-4 



f\Y\Z N ) 
/VIZ") 



(2.4) 
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Using (2.3) and the zero-one variable z /| ^=1 i A N ) 



where Z* represents N data sets, y of the data (and thus 55 
induces a partition of the data), T5 
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Minimize £ e,,...^,...^ 

Subject to £ 
*3 -'a 



tfi = t Mi) 



12.51 



P(r -ylZ") is the posterior probability of a partition y true 
given the data Z" The term partition is defined below; 65 
however* this framework is currently sufficiently general to 
cover set packings and coverings. 



iliy ...i N 
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-continued 

l M p and p = %...N-\) 

Zi x € 10, 1| for all ii i N . 

whexe 0 is arbitrarily defined to be zero. Here, each 
group of sums in the constraints represents the fact that 
each-non-dummy report occurs exactly once in a t4 track of 
data." One can modify this formulation to include multias- 
signments of one. some, or all the actual reports. The 
assignment problem (2.5) is changed accordingly. For 
example, if Zj * is to be assigned no more than, exactly, or no 
less than n (i * times, then the "=1" in the constraint (2.5) is 
changed to ^. =, ^ n,-*.** respectively. Modifications for 
group tracking and multiresolution features of the surveil- 
lance region will be addressed in future work. In making 
these changes, one must pay careful attention to the inde- 
pendence assumptions that need not be valid in many 
applications. 

Expressions for the likelihood ratios L,- t i*. can be 
found in the work of Poore. These expressions include the 
developments of Reid. Kurien. and Stein and Blackman. 
What's more, they are easily modified to include target 
features and to account for different sensor types. In track 
initiation, the N data sets all represent reports from N 
sensors, possibly all the same. For track maintenance, we 
use a sliding window of N data sets and one data set 
containing established tracks. The formulation is the same as 
above except that the dimension of the assignment problem 
is now N+l. 

Overview of the Lagrangian Relaxation Algorithm. 
(ABP) 

Having discussed the N-dimensional assignment problem 
(1.1), we now turn to a description of the Lagrangian 
relaxation algorithm. The algorithm will proceed iteratively 

for a loop k=l N-2. At the completion, there remains 

one two-dimensional assignment problem that provides the 
last step which yields an optimal (sometimes) or near- 
optimal solution to the original N-dimensional assignment 
problem In step k of this loop (summarized in Section 4) 
one starts with the following (N-k+l)-dimensional assign- 
ment problem with one change in notation. If k=l. then the 
index notation lj and L A are to be replaced by i x and M A , 
respectively. 
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j Subject to 

I = = ! U. 

Im-'n 

I = i-^-i 

for i p = 1 M p and p = Jt+ 2, ... , N - 1, 

15 I <i£ , -*- 1 -<-- 1 

To ensure that a feasible solution of (3.1) always exists for 
20 is a sparse problem, all variables with exactly one nonzero 
index (i.e.. variables of the form z ll0 0 *~** fori*. L A 

and V . . w . . . o""*" 1 for V= L ' * K ^ P= k+12 * • • ■ 
.N) assumed free to be assigned with the corresponding cost 

coefficients being well -defined This assumption is valid in 

25 the tracking environment. 

Subsection 3. 1 presents some of the properties associated 
with the Lagrangian relaxation of (3.1) based on relaxing the 
last (N-k)-sets of constraints to a two-dimensional one. A 
new approach to the problem of recovering a high quality 
feasible solution of the original (N-k+1) -dimensional prob- 

30 lem given a feasible solution (optimal or suboptimal) of the 
relaxed two-dimensional problem is described hereinafter. A 
summary of the relaxation algorithm is given hereinafter, 
following which we establish the maximization of the 
Lagrangian dual (an important aspect of the relaxation 

35 procedure) to be an unconstrained nonsmooth optimization 
problem and then present a method for computing the 
subgradients. 

The Lagrangian Relaxed Assignment Problem. The 
(N-k+l)-dimensional problem (3.1) has (N-k+1) sets of 
constraints. A (M r +l)-dimensional multiplier vector (Le., yf 
e R Mr 1) associated with the p-th constraint set will be 

denoted by u p =(u 0 p .u 1 p u M p r with u/sO being fixed 

for each p=k+2 N and included for notational 

convenience only. Now. the (N-k+l)dimensional assignment 
problem (3.1) is relaxed to a two-dimensional assignment 
4 * problem by incorporating (N-k-1) sets of constraints into 
the objective function via the Lagrangian. Although any 
(N-k-1) constraint sets can be relaxed, we choose the last 
(N_k_l) sets of constraints for convenience. The relaxed 
problem is 
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<l>*- 4 +i<« 1+2 u N )= Minimize <p N . l+ ii^' i * i ;u k * 2 ...u N )s 

Mini™*, £ «-*<&U ♦ 



(3.2) 



Minimize V ♦ £ £ 

*=i-2*,«0 

SubjecMo £ ^uVw-'*- 1 ^- 1 ' 



One of the major steps in the algorithm is the maximization 

of <&N-to.i(u M u") with respect to the multipliers 

(u* +2 u*). It turns out that is a concave, 

continuous, and piece wise affine function of the multipliers 

(u* +2 u*), so that the maximization of 4> yv _ Jt+A is a 

problem of nonsmooth optimization. Since many of these 
algorithms |???| require a function value and a subgradient 

of ^-jfc+i at any required multiplier value (u* +2 u*), we 

address this problem in the next subsection. We note, 
however, that there are other ways to maximize <$>n- M and 
the next subsection just addresses one such method. 
Properties Lagrangian Relaxed Assignment Problem 
For a function evaluation of O^+l, we show that an 
optimal (or suboptimal) solution of this relaxed problem 
(3.2) can be constructed from that of a two-dimensional 
assignment problem. Then, the nonsmooth characteristics of 
®N-**t arc addressed, followed by a method for computing 
the function value and a subgradient. 
Evaluation of <b N „ k+l . Define for each (l^u+i) an * ndex 

0* + 2- - • v jwHjM(UJ*.i). * • • JaXWi)) and a new cost 
function c. L 2 by 

org ixiin* 



Then. 

25 4. N . ttl ,U^ «)- 

Minimize (z 2 ; «* +2 u N ) 

''k 



(3.4) 
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Z Z fyk+A'M - Z Z w <' 

* t =0 



40 As an aside, two observations are in order. The first is that 
the search procedure needed for the computation of the 
relaxed cost coefficients in (3.3) is the most computationally 
intensive part of the entire relaxation algorithm. The second 
is that a feasible solution z N ~ k+i of a sparse problem (3.1) 

45 yields a feasible solution z 2 of (3.4) via the construction 



(3.3) 



7 1 



i. tf CL-'w = 1 for somc <'*♦«.». 4»k 

0, otlierwise. 



i,=0, 
1 M F and p « * + 2, 
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-4 



for </*,f 4+ i)*(0.0) 

<t+2-W W V 



thus, there are generally solutions other than the one nonzero 
index solution. 
The following Theorem 3. 1 gives a method for evaluating 
55 ^amh-i and states that an optimal solution of (3.2) can be 
computed from that of (3.4). Furthermore, if the solution of 
either of these two problems is e-optimal, then so is the 
other. The converse is contained in Theorem 3.2. 
Theorem 3.1. Let <o 2 dimensional assignment problem 
60 (3.4) and define to*"** 1 



(3.5) 



Given an index pair {\ k S M )Si M In) need not ^ 65 

unique, resulting in the potential generation of several 
subgradients (3.11). 



and fi+i>*«U 
if 0m,...I»)*Um J*) 
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-continued 

and (0,( 



t+1 r- 1 



*«t+2 



pmt+2 



..k+2 



Then ar^" 1 v problem and v 1 
<JWiyto. d* +2 \>*A addition, co, co^ 1 O^^yi/"^ i^A 
<^iY^ +2 ^* 



. . . .V 

„JV A 



Proof. Let to 



co 2 <J>*_*+iYCO J 



\TA <D 



YN-*+i Z °*" t ' 2 vA ° objective function values of (3.2) and (3.4). 
respectively. Direct verification shows that co +1 in (3.2) 
and ^al^TO)*^ 1 v k * 2 v" Q^+tfCif v k + 2 . . . .^remain- 
der of the proof, assume that co* y A 
Letx^YA 



i k _ 2 ...i N 

*oo = l if (k,k+i) = (0.0) and 

N 

Cm££...i N + 2 "O' -° for 501116 (i * +2 

/»»i+2 

=0 if (/ t ,it + i) = (0,0) and 

+ X < > 0 for 211 (i4+2 

^-*+2 



(3.6) 



ago = 1 if (/*, = (0, 0) and 

p-*+2 

u-Jo «= 0 if (4, « (0, 0) and 
ft 

+ £ uf p > 0 for all [i M i*>. 
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Note that x 2 y A .yx""** 1 ^ 2 v N A= <*Wi(r *>* +2 

r" - * 1 " 1 y A co co*"** 1 optimal solution of (3.2). Next 

o <l> A r- Jt+1 (£0 2 \)** 2 \T A co y A 

With the exception of one equality being converted to an 
inequality, the following theorem is a converse of Theorem 
3.1. 

Theorem 3.2. Let co""*" 1 y A and define to 2 
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The Nonsmooth Optimization Problem. Next, we address 
the nonsmooth properties of the function <J>A^ +1 yv)* += v"A 
explained in the following theorem, 

Theorem 3.4 Let y A lY*"*" 1 ) be the object 

function value of the (N-k+l)-dimensional assignment 
problem in equation (3.1) corresponding to a feasible solu- 
tion z"-* +1 of the constraints in (3.1). and let z^ +1 be an 
optimal solution of (3.1). Then. <J>am^iYU*" 2 is P iece " 
wise affine. concave and continuous in {v k * 2 v" and 



10 



Most of the algorithms for nonsmooth optimization are 
based on generalized gradients, called subgradients. given 
by the following definition for a concave function. 

Definition 3.5. At u=(u* +2 u*) the set 8<J> yv _^. 1 (u) is 

called a subdifferential of <t> s _ k ^ and is defined for the 
concave function O w _ A+1 (u) by 

i(u)£g T (Q)-u) for all coeR*"^ 1 x . . . jrR*** 

where g o p =ca o p =u o r =0 are all permanently fixed. (Recall that 
these were used for notational conveience only.) A vector 
gedO^ +1 is called a subgradient. 

There is a large literature on such problems, and we have 
tried a variety of methods including subgradient methods 
and bundle methods. Of these, we have determined that for 
a fixed number of nonsmooth iterations (e.g.. twenty), the 
bundle trust method Schramm and Zowe provides excellent 
quality solutions with the fewest number of function and 
subgradient evaluations, and is therefore our currently rec- 
ommended method. 

An Algorithm for Computing ^-k^-i ^ a Subgradient 
Most current software for maximizing the concave function 
^jv-a+i requires the value of the function and a subgradient 

at a point (u** 2 u*). Based on the previous two sections, 

this can be summarized as follows. 

1. Starting with problem (3.1). form the relaxed problem 
(3.2); 

2. To solve (3.2) optimally, defined the two dimensional 
assignment problem (3.4) via the transformation (3.3); 

3. Solve the two-dimensional problem (3.4) optimally; 

4. Reconstruct the optimal solution, say w* - *"** 1 . of (3.2) 
via equation (3.5) as in Theorem 3.1; 

5. Compute the function 



(3.10) 
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114 I 

p-i+2 /,-0 



x>"A co^ 1 y A co 2 (3.4), <IW,yto"-^ 
~»-*.i«o 2 v^ 2 v s A and (v^ 2 
V A 



Then to 2 y A O a 



v** 2 v N A= 
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M2 



6. A subgradient is given by 



(3.11) 



and 



Proof: Let co^ 1 co 2 <J>^j(co ~ ^ 

O^^co 2 v** 2 v* r A objective function values of (3.2) and 60 ^) m 



- /v-*-hiv~ - - — - -j — v 

(3.4). respectively. Direct verification shows that © (3.4) 
and ^^^yu/^** 1 \)*^ 2 d^A^^^^co 2 d** 2 d n A remain 



(0 2 (0 



co 



and ^^yoT-* 1 - 1 ^AlO^^co 3 iTA re: 
der of the proof, assume that or^ y A and construct 
cu""* 1 " 1 _ by replacing or*"** 1 y A co""** 1 - 
<t>N-*+ifN-**lv v N A=0^ 1 ((o 2 v k + 2 v"A* 
QM^ifsf 1 *** 1 '* v** 2 v N A (jf-** 1 inequality is in fact an 
equality. This proves the theorem. 
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for i p = l ? ...,Mp and p » 2, ... , 

Several remarks are in order. First, the optimal solution of 
the nunimiation problem (3.2) is required before one can 
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remove the minimum sign, replace z^ +i by the minimizer 
w *-a-*i and differentiate with respect to u, r to obtain a 
subgradient as in (3.11). If w AMc ' 1 ' 1 is an approximate solu- 
tion of (3.2), then the subgradient and function values are 
only approximate with accruacy depending on that of 
w* - ***" 1 . Although one can evaluate the sums in (3.10) and 
(3.11) in a straightforward manner, another method is based 
on the following observation. Given the multiplier vector 
u,„ 2 u„), let {(Wl^^a^id^i)}^^ 1 ^. be an enu- 
meration of indices ^lA+i} of w * ( or me 2 ^ ccs of 
w"~* +1 constructed in equation 3.5)) such that 
^kuww- 1 for (MWUUMOfl) and (1,(1^). 
ijfc+i(U +l ))=(0,0) for 1 A+1 =0 regardless of whether Wqo =1 or 
not. (The latter can only improve the recovered feasible 
solution.) The evaluation of the bracketed quantity in (3.11) 

for a specific i,,>=l and p=k+2 N is one minus the 

number of times the index value i p appears the (p-k+1)"' 
position of the (N-k+l)tuple in the list ({1^+1)^+1(^+1)' 

Finally, we nave presented a method for computing one 
subgradient. If w vv ~* +1 is the unique optimal solution of 
(3.2), then 0 Af _ A+1 (u) is differentiable at u, g is just the 
gradient at u. and the subdifferentiai o4>yv_ t+1 (u)={g} is a 
singleton. If, corresponding to u.w*"*'** 1 is not unique, then 
there are finitely many such solutions, say 
{w^***(l), . . . ,w JV ~** +A (m)} with respective subgradients 

jg(l) g(m)|, the subdifferential50(u)is the convex hull 

of {g(l) g(m)}. These nonunique solutions arise in two 

ways. First, if the optimal solution of the two dimensional 
assignment problem is not unique, then one can general all 
optimal solutions of the two dimensional assignment prob- 
lem (3.4). Corresponding to each of these solutions and to 
each index pair (1*4^)01 each solution, compute the indices 

(k +2 ,j^) in (3.3). If these j*s are not unique, then we can 

enumerate all the possible optimal solutions $/ N ~ k + l of (3.2). 
Given these solutions, one can then compute the correspond- 
ing subgradients from (3.11). 

In most nonsmooth optimization algorithms, one uses an 
epsilon-subdifferential 

Definition 3.x. At u=(u* +2 u") the set 8 c <I> A ^ k+1 (u) is 

called an epsilon subdifferentiai of <IV-*+i and * s defined for 
the concave funtion Q>N-k+i( u ) 

i(u)Sg T (v>- U y* for all wR^'x . . . xR K * +l }. (3.9) 

where g o p =GV=u o /, =0 are all permanently fixed. (Recall that 
these were used for notational convenicence only.) A vector 
g€t)4V_ A+1 (u) is called a subgradient. 

HOW DO YOU COMPUTE gea<tWi(u) 
The Recovery Procedure. The next objective is to explain 
a recovery procedure, i.e.. given a feasible (optimal or 
suboptimal) solution co 2 7 A 7 <jf~ k+L 7 A via Theorem 3.1), 
generate a feasible solution z""*** 1 of equation (3.1) which 
is close to (o 2 specified. We first assume that no variables in 
(3.1) are presassigned to zero; this assumption will removed 
shortly. The difficulty with the solution to""** 1 <o it need not 
satisfy the last (N-k-1) sets of constraints in (3.1). (Note, 
however, that if (O 2 (3.4) and (0 A/ ~* +l (o relaxed constraints, 
then ay*"** 1 7 A A o) recovery proceudre described here is 
designed to preserve the 0-1 character of the solution u) 2 7 60 
A If to Vl+| 2 =l and 1**0 or 1^*0, the corresponding feasible 
solution + z *-* +l of (3.1) is construed so that z v . . y> . . . 
o N-A+1 =l for some (i M , . . . i w ). By this reasoning, variables 
of the form z^^ iJ*~** X can be assigned a value of one 
in the recovery problem only if C0oo 2 =l. However, variables 
AxMt, i N ~ k * 1 wiU 1* treated differently in the recovery 
procedure in that they can be assigned either zero or one 



independent of the value of (Ooq 2 . This increases the feasible 
set of the recovery problem, leading to a potentially better 
solution. 



Let {(lA^^id^)} 
indices {1 



CD, 



i A+l (i A+1 )MO,0) for 1* +1 =0 regardless of whether (Ooo*= 1 or 
not. To define the N-k dimensional assignment problem that 
resores feasibility, let 



£, L +l } of or 7 or~* +i equation (5)) such 
*=1 for (Wl fcfl ),Ul w )M0,0) and (1^(1 



be an enumeration of 
that 



*n = «fi 



for k = 1 



(3.12) 



" c '*Wt + iWi + 2- *n 

" C 'l«'2...*+i> i 2*'2...1+i^ l '»...i+l> - 
'*('u+l>'t+l(4t+lV"l+2-" l N 

for **2 and U = 0, 



where 



• oW M Y**JUi( • ■ ■ ('*('*♦!» ■ • • )X3-13) 
for m=2 k+1 and where o denotes function composition. 

For example ^^(U) and ^234^2° WU^WWU))* 

Then the (N-k)-dimensional assignment problem that 
restores feasibility is 



Minimize 



Z 

'i+l'i+2" 



r N-k 



N-k 



in 



(3.14) 



Subject to Yj 



1, Im = 1 *4+l. 



z 



• I, *i + 2 = I. 



, M* + 2. 
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Z 



*N 



''i+l'i+l-'W 



= 1, 



for i p = 1 M p and p = K + 3 N-L 
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'i+i'i+i—'w-i 



Let Y be an optimal or feasible solution to this N-k 
dimensional assignment problem. The recovered feasible 
solution z N is defined by 



(3.14) 



50 



1 if h « Mfc... a + i). h - i2('2..a + i>. h = fj('2-..i + l). 



otherwise. 



1 <Vfi('i+i> and *i 1+ ,/ t+2 - 'w 



= 1. 



55 



Said in a different way, the recovered reasible solution z* 

corresponding to the multiplier set { u* +2 u N ] is then 

defined by 



E 'H'2-..4+l^<'2...i+|V3tfj...i + |Viatt + |Vi + i(/A + |Wi + j-^ ' 



1 if y, 



Wi+2" 'W ' 
otherwise 



65 



where l m . . . 4*1 is defined in (3.13) and o denotes function 
composition. 
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The next objective is to remove the asumption that all cost 
coefficients are defined and all zero-one variables are free to 
be assigned. We first note that the above construction of a 
reduced dimension assignment problem (3.14) is valid as 
long as all cost coefficients c*~* are defined and all zero-one 
variables in z*~* are free to be assigned. Modifications are 
necessary for sparse problems. If the cost coefficient c w/fcti) 

A 



N ~ jw ' 1 is well-defined and the zero-one variable 



m n v i s free to be assigned to zero or one, 10 

then define c w - +2 „"-*=c, . « J*** 1 as in 



i"-* being free to be assigned. 



(3.12) with z w „ . 

Otherwise. z Ww . . . ^ is preassigned to zero or the 
corresponding arc is not allowed in the feasible set of arcs. 
To ensure that a feasible solution exists, we now only need 15 
ensure that the variables z, i+i0 , . . o""* are free to be assigned 
for 1 A+1 =0.1. . . . with the corresponding cost coeffi- 
cients being well-defined. (Recall that all variables of the 
form z tk0 . . . o"-** 1 and Zo . . . 0ip o . . . o""** 1 ™ f £^(to be 20 
assigned) and all coefficients of the form c li0 0 
Co oy> o^-*" 1 « well-defined for 1*=0. . . . X* and y= 
0. ... M p for p=k+l. . . . JSL) This is accomplished as 
follows: If the cost coefficient c lk(l y ( i M yo . . . o*""** 1 ^ 
well-defined and 2, (WW/i+l)0 of-**' is free, then define * 
c iM o . . . o^=WwW>o - • • o^ 1 with z, +l0 ... being 
free. Otherwise, since all variables of the form z ik0 0 
and Zo, o o*"*"*" 1 are free to be assigned with the 
corresponding costs being well-defined, set c 1m0 0 = 
Wo ■ • ■ o^+cowwo • ■ - o^ 1 where the firstterm is 
omitted if 1*(W=0 and the second, if lw^iH- For 
WwH and i^Ot+i)^ define c 0 . . . 0 *-*=Co . . . . 

The last step k=N-2. The description of the algorithm ^ 
ends with k=N+2. The resulting recovery problem defined in 
(3.14) with k=N-2 is 



1 if IVi'V^ 
0 otherwise. 



where l m ^ is defined in (3.13) and o denotes function 
composition. To complete the description, let {(ln-iflv)^ 
(Wh ^0 °e an enumeration of indices {l/v-ii^} of Y such 
that ^ j(WJj ^f»l for and (WW>. 

i^WHO.O) for 1,^=0 regardless of whether Y^l or not 
Then the recovered feasible solution can be written as 



1 if t] = hik-N-lX '2 = '2^2- Af- A *3 = *j('s-*-l). 

0 otherwise. 
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The Upper and Lower Bound The upper bound on the 
feasible solution is given by 



V*(£")=I ij . . . i/f C *! - - . V **1 ■ ■ • V 



(3.21) 



and the lower by dvy 3 . . . or*) for any multiplier value 
(u 3 u*). In particular, we have 



O^ir* iOS^ OS V^f) 



(3.22) 



<h s Minimize £ £ ^-iijv^-i'w s V2(=2) 
Subject to ^ ^N~liN'^ i am -J— iw-i ' 

^-iV € 10,11 foraU 



where (u 3 u*) is any multiplier value. C 3 . . . r*A 

niaximizer of ^y 3 .... u*). ~" multidimensional assign- 
ment problem (3.1) and l N is the recovered feasible solution 
defined by (3.20). Finally, we conclude with the observation 
that Vjv(2)=V 2 (Y) where Y is the optimal solution of (3.17) 
used in the construction of z in equations (3.18H3.20). 

Reuse of Multipliers. Since the most computationally 
expensive part of the algorithm occurs in the maximization 
of 
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Let Y be an optimal or feasible solution to this 
2-dimensional assignment problem- The recovered feasible 
solution z* is defined by 



starts or multipliers close to the optimal are fundamentally 
important for real-time speed. The purpose of this section is 
to demonstrate that the multiplier set obtained at stage k£ 1 
provide good starting values for those obtained at step k+1. 

Theorem 3.xx. Let (u 3 u") denote a multiplier set 

obtained in the maximum of O^yy 3 . . . .u*). Then this 
multiplier set satisfies the string of inequalities 

<D^ ^SO^u 4 , . . - */*)£ . . - S<I>,(OS0 2 S 

where after the first step in the niaximization of 4> N multi- 
pliers are not changed in the remaining steps. Furthermore, 



(3.18) 



55 



1 if ii « <i«2» w-i)» h - hih~N-\\ hih -ft-i) 

*W-2 ° f*-2('Af-2JV-lV'W-lt'w-l) ^ 1 

or if (I*-i.4v)°<0,0) 
0 otherwise. 



60 



to improve this inequiabty, let (u 
denote a maximizer of Oy^^u** 2 .u"). Then, we have 

& * Wi-^** 1 ^ t3 - ra) 



Inequalities depend on solution of 2d assignment prob- 
65 lem* 

Said in a different way, the recovered feasible solution is Proof: Suppose we have a value of the multipliers 
then defined by (u** 2 Ji") obtained in the maximization of 
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*w-i + i(«* +J . -M «*> s Minimize «** 2 «*) 

Subject to: £ ^^..^-ia-i...^. 

4t+l'i+2"'W 
Vi+2"'* 



where 



78 

-continued 

'i+2 1, .... A/i+j. 



Since the constraint I, 



10 



iaT-* +, =1 for 1^2=1.. 
so that one has 



.Mjk2 is now imposed, the feasible region is smaller. 



EE 



i/^^.O^ 3 «* 

where the fact that Q^+i does not depend on the multiplier 
i 5 set u* +2 due the added constraint set. Also, the last equiva- 
lence follows from the fact that ^^-nY 4 * 3 • • • .u N ) is the 
relaxed problem (3.xx) (with k replaced by k+1) for the 
recovery problem (3.xx) in step k of the above algorithm. 

Continuing this way. one can compute (u 3 u") at the first 

20 step (k=l), fix them thereafter, and perform no optimization 
at the subsequent steps to obtain 



These need not be the maximizer; however, we do assume 
that we have solved the above minimization problem opti- 25 

mally to evaluate <I>/vr_* + i(u* +2 u"). Just as in the 

definition of the earlier recovery problem, let {(i k (i k+i )d k+l 
(U+i)}/^/** 1 ^ be an enumeration of indices {\ k S k+l ) of the 
optimal solution w 2 of problem (3.4) (or the first 2 indices 
of the solution w^* 1 of the relaxed problem (3.2) con- 30 
structed in equation (3.5)) such that w, i(/i t)ii l(/il ) 2 =l for 
(WUWUW0.0) and (lA.i^ifUW.O) for 
ljfc +1 -0 regardless of whether Woo 2 =l or not. 

Then, the final value of the objective function can be 35 
expressed as 

<JWi(«** a «")sMinimizc ^-^(z^ 1 ;^ 2 , ...,«*) 

Subject to t lM . . . tM z t ^ lVkfli i MVM . . . ^ . . . xA+i""** 1 * < 3m ) 40 

where 



<t>^u 3 «")£<IW« 4 . . • 3* 3 («*)S<MV' W (2) 

where <J> 2 Y & 

To explain how to improve this inequality, let 
( u 2 + *jv-a+i yA/jv-^ij denote a maximizer of <tV_* +1 

(u* +2 u"). Then by the same reasoning one has (3.xx) 

as stated in the theorem. Q.E.D. 
Summary of the Lagrangian Relaxation Problem. (APB) 
Given the multidimensional assignment problem 



Minimize £ 'V 'jv^i 



(4.1) 



(3. xxx) 



z 



'*+l'i+2-'* 



P«t+2 



If now another constrain set, say the (k+2) ,A set, is added to 
this problem, one has 



45 



Subject to ^ Zi r . .i N = 1 (ij = 1 M2) 

Z 

(l p « 1 Mp and /> = 2, ... , A/- l) 

Z = 1 B 1 m *> 

'l*2 "'W-l 

-ijycf0.il fee «U 



<*W+iU*** 2 . ..-«") a Minimize $N-k*\K*?~ k * X \ « i+ 



Minimize 



z 



-Nil 



<* Z 
p-t+2 



Subject to £ 

<i+2"'* 



i^ + ifi+i)"» Z Z"£ 



,1.4*1 1 U+l. 



where ^ /is arbitrarily defined to be zero and is included 

50 for notation al convenience and where the superscript N on 
both c and z is not an exponent, but denotes the dimension 
of the subscripts and the assignment problem as stated in 
(2.5). In relaxed and recovery problems c 0 0 N need not be 
zero! In this problem, we assume that all zero-one variables 

55 z ( , * with precisely one nonzero index are free to be 
assigned and that the corresponding cost coefficients are 
well-defined. (This is a valid assumption in the tracking 
environment.) Although not required, these cost coefficients 
with exactly one nonzero index can be translated to zero by 

60 cost shifting without changing the optimal assignment. 
Having explained many of the relaxation features, it is 
now appropriate to present the Lagrangian relaxation 
algorithm, which iteratively constructs a feasible solution to 
the N-dimensional assignment problem (4.1). 

65 Algorithm 4.1 (Lagrangian Relaxation Algorithm). To 
construct a high quality feasible solution, enoted by of 
the assignment problem (4.1), proceed as follows: 
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0. Initialize the multipliers (u* +2 u"). e.g.. (u 

irXO.' . • • .0). & ft £ (51) 

Minimize 2, 2- L' 1 * 1 ***** 

For k- 1 N-2. do < , -o #2 -o < 3 -0 

1. For the Lagrangian relaxed problem (32) from the 5 "2 "2 

problem (3.1) by relaxing the last (N-k-1) sets of con- Srf » ecI to fafe* 1 ™ = *' " = " 

2. Use a nonsmooth optimization technique to solve £ 2 ■= h - 1 ^ 

Maximize {O^^jhycR**^ 1 for p=k+2, . . . Jff with V=0 being w , «, 

feed} < 4 - 2 > £ £ = U 3 = 1 "3 

I j eO^oO 

where ^-^(u^ 2 u*) is defined by equation (3.2). The e {0, n for all i,; 2 ij. 

algorithm following problem (3.9) provides one method for 15 

computing a function value and a subgradient out of the m t r ... . c /e iV . • ♦ * 

^ 0 , , v • j - To ensure that a feasible solution of (5.1) always exists for 

subdifferential at (u - u ). as required in most a spme ^ variable s with exactly one nonzero 

nonsmooth optimization techniques. index variables of the formz^ .z^. and z^ for 

3. Given an approximate or optimal maximizer of (4.2). 2Q 1 M p and p=123 are assumed&ee to be assigned with 

say (G** 2 0"). let w 2 assignment problem (3.4) the corresponding cost coefficients being well-defined. This 

corresponding to this maximizer of fl)^ 3 u"). assumption is valid in the tracking environment [**refPoa. 

* **refPob]. 

4. Formulate the recovery (N-k)-dimensional problem Lagrangian Relaxed Assignment Problem. The three 
(3.13). modified as discussed at the end of subsection (3.3) dimensional assignment problem (5.1) has three sets of 
for sparse problems. At this stage. l N as defined in (3.15) 25 constraints and one can describe the relaxation by relaxing 

contains the alignment of the indices i*+i}- any of the three sets of constraints, the description here is 

Formulate the final two-dimensional problem (3.17) and based on relaxing the last set of constraints A <M 3 +1)- 

complete the final recovered solution z" as in (3.20). dimensional multiplier vector (i.e.. u 3 eR ) associated 

. 30 with the 3 th constraint set will be denoted by u~=(u 0 . 

To explain the lower and upper bounds, let <J>„ defined in JU 3 u M Y with u 0 3 sK) being fixed for notational 

(3.2) with k=l. let V^z") be the objective function value of convenience. (The zero multiplier u 0 3 is used for notational 

the N-dimensional assignment problem in equation (2.5) convenience.) Now. the three dimensional assignment prob- 

correspondingto a feasible solution z* of the constraints in 1^ (5 j s re laxed to a two-dimensional assignment prob- 

(2.5). and let z 7 *""** 1 (2.5). Then. 3$ iem by incorporating last set of constraints into the objective 

function via the Lagrangian. Although any constraint set can 

&^u* « w )gv A <?')gvv^) (43) be relaxed, we choose the last set of constraints for conve- 
nience. The relaxed problem is 

is the desired inequality. 



40 



M] M 2 M; (5.2) 

3 <tf*) = Minimize <fo{f ;u*)e Y Y ^4*i''3^'J + 
COMMENTS n -0 i 2 -0 * 3 -0 



1. Step 2 is the computational part of the algorithm. 
Evaluating .O^^j search procedure requires 99% of the 45 { - m0 
computing time in the algorithm. This part uses two dimen- 
sional assignment algorithms, a search over a large number & Minimize Y Y Y (c 3 + ]z? - Y "? 
of indices, and a nonsmooth optimization algorithm It is the B fciip>£& Wz 13 ^ & ' 
second part (the search) that consumes 99% of the compu- 
tational time and this is almost entirely parallelizable. 50 subject to « 1 ,• + 1 = 1 Mi 

2. In track maintenance, the warm starts for the initial 
multipliers for the first step are available. For the relaxation 



; 2 -Oi 3 -o 

procedure, initial multipliers are available in step 2 from the 
prior step in the algorithm, 55 



: 1, I + 1 « 1 M 2 



3. There are many variations on the above algorithm. For 0 ne of the major steps in the algorithm is the maximiza- 

example, one can compute a good solution on the first pass y on c f <j> 3 t 3 A u 3 . It turns out that <D 3 piecewise affine 

(k=l) and not perform the optimization in (2) thereafter. This function of the multiplier vector u 3 , so that the maximization 

yields a great solution. Thus one can continue the optiroi- 60 of <I> 3 Since many of these algorithms require a function 

zation at the first pass, and immediately compute quality value and a subgradient of <D 3 y 3 A address this problem in 

feasible solutions to the problem. the next subsection. We note, however, that there are other 

ways to maximize <Z> 3 subsection addresses but one such 

Lagrangian Relaxation Algorithm for the 3-Dimensional method. 

Algorithm. (ABP) In this section, we illustrate the relaxation 65 Properties Lagrangian Relaxed Assignment Problem 

algorithm for the three dimensional assignment problem. For a function evaluation of 0 3 y suboptiroal) solution of 

Having discussed the general relaxation scheme, this relaxed problem (5.2) can be constructed from that of a 
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two-dimensional assignment problem. Then, the nonsmooth 
characteristics of <t> 3 followed by a method for computing a 
subgradient. 

Evaluation of <t> 3 . i 2 ) j 3 =j3(ii.i 2 ) and a new cost 
function c (i ^ 2 by: 5 

h s A('i. k) = arg mini*?, V3 + «? 3 I is - 0, 1 My \ < 53 > 

d} ii2 =c5 3 for (ii./2>*<0.0) 



*3 



"*t'i = Z *Wi for (/| ' ,j) * (0 * 0) 



(5.7) 



« 3 



10 



cJo = Z Minimum (0, cj^ + «J 3 | 

13 c0 



Given an index pair (i^)* j 3 need not be unique, resulting 15 
in the potential generation of several subgradients (5.11). 
(We further discuss this issue at the end of Subsection 5.2.3 .) 
Now, define 



4*m 3 ) = Minimize ^Mz 2 ; " 3) s Z Z ^ib^i'i " Z 



(5.4) 



20 



Subject to £ jj |l2 = 1, 
i 2 =o 



25 



woo 2 =0 if (i.,i 2 >=(00) and c^+u^ for all i 3 
Woo 2 =l if (i,,i 2 >=<00) and coq^+u^^O for some i 3 . 

Then w 2 is a feasible solution of the problem (5.4) and 
<J> 3 (w 3 ;u 3 )^<t> 3 (w 2 ;u 3 ). If, in addition, w 3 is optimal for 
(5.2), then w 2 is an optimal solution of (5.4), <l> 3 (w 3 ;u 3 )= 
<J) 3 (w 2 ;u 3 ) and <J> 3 (u 3 )=0 1 (u 3 ). 
The Nonsmooth Optimization Problem 
An Algorithm for Computing <1> 3 and a Subgradient 
Given u 3 the problem of computing <J> 3 y 3 A subgradient of 
4> 3 y 3 A 

1. Starting with problem (5.1). form the relaxed problem 
(5.2) 

2. To solve (5.2) optimally, define the two dimensional 
assignment problem (5.4) via the transformation (5.3) 

3. Solve the two-dimensional problem (5.4) optimally. 

4. Reconstruct the optimal solution, say w 3 e y A equation 
(5.6) as in Theorem 5.1. 

5. Then 



Z^ = l > 1 

z 2 ^ e(0, 1| for all ij, / 2 - 



*#( w 2 M 3 



(5.10) 



30 



As an aside, two observations are in order. The first is that 
the search procedure needed for the computation of the 
relaxed cost coefficients in (5.3) is the most computationally 
intensive part of the entire relaxation algorithm. The second 35 
is that a feasible solution z 3 of a sparse problem (5.1) yields 
a feasible solution z 2 of (5.4) via the construction 



"3 
i 3 =0 



Z Z *'1'2'3 ■ 



6. A subgradient is given by substituting w 3 objective 
function (5.2), erasing the minimum sign, and taking the 
partial with respect to u ( 3 The result is 



l.fi 



if z?,^ = 1 for some (/,, i 2 , 13); 



40 



0, otherwise. 



4 = 



Z Z ^ " 1 



(5.11) 



for *3 = 



, A/3. 



Thus the two-dimensional assignment problem (5.4) has 
feasible solutions other than those with exactly one nonzero 45 
index if the original problem (5.1) does. 

The following Theorem 5.1 states that an optimal solution 
of (5.2) can be computed from that of (5.4). Furthermore, if 
the solution of either of these two problems is e e then so is 
the other. The converse is contained in Theorem 5.2. 5Q 

Theorem 5.1. Let co 2 dimensional assignment problem 
(5.4) and define w 3 by 

w^*^^ 2 if i,=j 3 and (i»,i 2 )*<0,0); 

"1^=0 ^ »3*b and (>i,i 2 M0,0); 55 

Wom 3 3 =1 »f <W+«ij 3 £<>; 

Wooi/^OlifW-Ki^. (5.5) 

Then w 3 is a feasible solution of the Lagrangian relaxed 
problem and O f 3 A O f 3 A e e 2 optimal for the two 60 
dimensional problem, then w 3 is an optimal solution of the 
relaxed problem and <I> 3 (u 3 )=<I> 3 (u 3 ). 

With the exception of one equality being converted to an 
inequality, the following theorem is a converse of Theorem 
5.1. 65 

Theorem 5.2. Let w 2 be a feasible solution to problem 
(5.2) and define w 2 by 



The Recovery Procedure. The next objective is to explain 
a recovery procedure, i.e., given a feasible (optimal or 
suboptimal) solution w 2 of (5.4) (or w 3 of (5.2) constructed 
in Theorem 5.1), generate a feasible solution z 3 of equation 
(5.1) which is close to w 2 in a sense to be specified. We first 
assume that no variables in (5. 1) are preassigned to zero; 
this assumption will be removed shortly. The difficulty with 
the solution w 3 in Theorem 5.1 is that it need not satisfy the 
third set of constraints in (5.1). (Note, however, that if w 2 is 
an optimal solution for (5.4) and w 3 , as constructed in 
Theorem 5.1, satisfies the relaxed constraints, then w 3 is 
optimal for (5.1).) The recovery procedure described here is 
designed to preserve the 0-1 character of the solution w 2 of 
(5.4) as far as possible: If w (|ii 2 =l and ij*0 or i 2 *0, the 
corresponding feasible solution z 3 of (5.1) is constructed so 
that (, 3 =1 for some (ii.i 2 J3> • By this reasoning, variables 
of the 'form z^ 3 , can be assigned a value of one in the 
recovery problem only if w^^l. However, variables z^ 3 
will be treated differently in the recovery procedure in that 
they can be assigned either zero or one independent of the 
value Woo 2 . TW S increases the feasible set of the recovery 
problem, leading to a potentially better solution. 

Let {(ii(l2)i 2 ( 1 ?.)}f I -0 ^ an enumeration of indices 
{ii4 2 } °f w (° r me m * st ^ indices of w 3 constructed in 
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equation (5)) such that Wi,^).*^)"— 

1 for (i^.i.^MO.O) 
and (ii(l 2 )i 2 (I 2 ))=(0.0) for 1 2 =0 regardless of whether 
Wqo^I or not. (The latter can only improve the quality of the 
feasible solution.) 

Next to define the two-dimensional assignment problem 
that restores feasibility, let 



Let Y be an optimal or feasible solution to this 
2-dimensional assignment problem. The recovered feasible 
solution z 3 is defined by 



i: 



if «i c iiihX h » »2(':> and ^'3 * 1; 



(5.15) 



0 otherwise. 



(5.12) 



Then the two-dimensional assignment problem that restores 
feasibility is 



10 



h "3 

Minimize £ £ 4*4'3 
^ -0/3-0 

Subject to ^ 4 13 = 1, = 1 La. 

<3-0 
in 

X 4o = x ' i3 = 1 M ^ 

4*3 €iC II for all / 2 ,/ 3 . 



(5.14) 



Said in another way. let {(UU^i^)}^*^ be an enumera- 
tion of indices {l 2 i 3 } of Y such that Y^ Vj</j) =l for (1 2 (1 3 X 
i 3 (l 3 )*(0.0) and (l 2 (l 3 )i 3 (l 3 M0.0) for* 1 3 =0 regardless of 
whether Y^l or not Then the recovered feasible solution 
can be written as 



15 



1 if '1 = iiOnl 12 = iiikiX ii = I3C3) 
0 otherwise 



(5.16) 



20 The upper bound on the feasible solution is given by 

«i JU 2 M 3 (5.17) 



25 

The next objective is to remove the assumption that all 
cost coefficients are defined and all zero-one variables are 
free to be assigned. We first note that the above construction 
of a reduced dimension assignment problem (5.13) is valid 
as long as all cost coefficients c : are defined and all zero-one 30 
variables in z 2 are free to be assigned. Modifications are 
necessary for sparse problems. If the cost coefficient 
^lihMW 3 * s we U-defined and me zero-one variable 



35 



and the lower by 0 3 y 3 A y 3 A 
In particular, we have 

where u 3 is any multiplier value, (u A 3 0 3 y 3 Ae z 3 assign- 
ment problem and z 3 is the recovered feasible solution 
defined by (5.16). 

Other Relaxations for the Multidimensional Assignment 
Problem.(ABP) In this section, we briefly describe other 
possible relaxations and their implications. These include 
linear programming relaxations and the corresponding 
duals. 

Recall that one starts with the definition of the zero-one 



^ith&chv* IS t0 assigned to zero or one. then define 

'L^iMMUh 35 m < 5 " 1 ? with 2 *,ft>*W 3 bein S free t0 be 
assigned. Otherwise, z ¥j is preassigned to zero or the 

corresponding arc is not allowed in the feasible set of arcs. 

To ensure that a feasible solution exists, we now only need 

ensure that the variables z h0 2 are free to be assigned for vari ablez~" 7=1 and cost coefficients to form foe problem 

1 2 =0.1 with the corresponding cost coefficients being 40 @ v ' N 

well-defined. (Recall that all variables of the form z il00 3 . 
Zq^ 3 and Zq^ 3 are free (to be assigned) and all coefficients 
of the corresponding form c^^.c^q 3 and Cqq^ 3 .) This is 
accomplished as follows: If the cost coefficient c^^^ 3 is 



where Cq 0 is arbitrarily defined to be zero. Here, each 
group of sums in the constraints represents the fact that each 
non-dummy report occurs exactly once in a 4i track of data**. 
3 ...y-.. One can modify this formulation to include multiassign- 

well-defined and 2 y ^>o^s ^^*^^Jf 3° ^J^H 45 ments of one, some, or all the actual reports. The assignment 

" °* problem (2.5) is changed accordingly. For example, if z ik k is 
to be assigned no more than, exactly, or no less than n, 4 * 
times, then the "=1** in the constraint (2.5) is changed to 



(*2X> 3 w * m z t& being free. Otherwise, since all variables 



the form z l|0 3 and z^o are free to be assigned with the 



corresponding costs being well-defined, set c^ =c Ji(Ji)00 + 
Cq^jq 3 where the first term is omitted if ij(l 2 )=0 and the 
second, if i 2 (l 2 >=0- For ii(l 2 )=0 and i 2 (l 2 )=0. define Coo 2 = 
c ooo • 

The final recovery problem. The recovery problem for the 
case N=3 is 



Minimize £21^^* 



Subject to 2 4fj ™ ' 2 * *• *" * 



(5.14) 55 



60 



"S-ss.^nf */* respectively. Modifications for group tracking 
50 and multiresolution features of the surveillance region will 
be addressed in future work. In making these changes, one 
must pay careful attention to the independence assumptions 
that need not be valid in many applications. 

Again, the recent work of Poore and Drummond signifi- 
cantly extends the formulation of the track maintenance and 
initiation to new approaches. The discussions of this section 
apply equally to those formulations. 

The Linear Programming Relaxation. In the linear pro- 
gramming relaxation, one replaces the zero-one variables 
constraints 



4^ €{0 t l| for all fc.ij. 



z*i . . .j/10,1} for all i„ ...4* 
with the constraint 

OSz^ . . . ijv Sl for all i,, . . . i N 



(6.2) 



(6.2) 



65 

Then, the problem (6.1) can be formulated as a linear 
programming problem with the constraint (6.2) in (6.1) 
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replaced by (63) with a special block structure to which the 
Dantzing- Wolfe decomposition is applicable. Of course, 
after solving this problem, one must now recover the zero- 
one character of the variables in (6.1) and there are many 
ways to do this, such as using the two dimensional assign- 
ment problems. Commercial software is also available. 

The Linear Programming Dual and Partial Lagrangian 
Relaxations. Given the linear programming relaxation, one 
can formulate the dual problem or the partial Lagrangian 
relaxation duals with respect to any number of constraints. 
In particular, this is precisely what is done in Section 3 on 
the Lagrangian relaxation algorithm presented. The much 
broader class of algorithms provided in the U.S. Pat. No. 
5.537,119 (Aubrey B. Poore, Jr., Method and System for 
Tracking Multiple Regional Objects by MultiDimensional 
Relaxation. U.S. Pat. No. 5,537.119, filed Mar. 14, 1995. 
issue date of Jul. 16, 1996) can be modified to remove the 
zero-one character when one relaxes M sets of constraints to 
an N-M dimensional problem and recovers with an N-M+l 
dimensional problem. This avoids the problems associated 
with the NP-hard N-M and N-M+l dimensional problems. 
However, to restore the zero-one character, one can do it 
sequentially with an assignment problem or with one of the 
many zero-one rounding techniques. These formulations are 
easy to work out and thus the details are omitted. 

The complete dual problem is another way of solving the 
LP problem and may indeed be more efficient in certain 
applications. In addition, the solution of this dual problem 
may provide excellent initial approximations to the multi- 
pliers for Lagrangian relaxations. 

Multisensor Resolutions Problems and Their Solution via 
Linear Nonlinear Assignment Problems and Inequality Con- 
straints. (ABP) The mathematical development for inequal- 
ity constraints is understood by Aubrey B. Poore, but this 
work has not been written up. This feature is important in 
combining information from sensors with different resolu- 
tions such an IR and Radar. 

Formulation and Algorithms for the Distributed Sensor 
Tracking: Track-to-Track Correlations. (ABP) 

Hot Starts for Track Maintenance and Initiation: Bundle 
Modifications (ABP) 

Thus reuse of multipliers and the first proof that this reuse 
is actually valid was presented in this section. The reuse in 
track maintenance is demonstrated and discussed in the 
work of Poore and Drummond [xx]. The only item left is the 
modification of the bundle of subgradients for the use with 
the multiplier values as one goes through the recovery 
problem as in Section 3.6 or as one moves the window in 
track maintenance. It is this aspect of the non smooth opti- 
mization that adds an order of magnitude to the speed of the 
algorithms. This work is based on both a mathematical proof 
as in Section 3.6 as well as computational experience and 
heuristics. 

Adaptive Window Size Selection (ABP) 

These data association algorithms, based on the multidi- 
mensional assignment problem, range from sequential pro- 
cessing to many frames of data processed all at once! The 
data association problem for multisensor multitarget track- 
ing is formulated using a window sliding over the frames of 
data from the sensors. This discussion focuses on. the work 
of Poore and Drummond and not on the formulations in 
Section 2. Firm data association decisions for existing track 
are made at, say frame k, with the most recent frame being 
k+M. Decisions after frame k are soft decisions. Reports not 
in the confirmed tracks are used to initiate tracks over frames 
numbered k-I to k+M. 
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The length of these windows varies from sequential 
processing, which corresponds precisely to 1=0 and M=l, to 
user defined large values of I and M. (In the case of 
sequential processing, we also have a temporary track file of 

5 dynamically feasible tracks, but incorrect data association.) 
There is a marked change in performance over this range. As 
the two windows become exceedingly long, there is little 
statistically significant information gained and indeed per- 
formance can degenerate slightly. Thus, one must manually 

10 find the correct window lengths for performance in a given 
scenario and the algorithms do not change for a changing 
environment. Thus we propose an adaptive method for 
adjusting the window lengths. (The method has been highly 
successful for selecting the correct window length for mul- 

15 tistep methods in ordinary differential equations.) 

1. Compute the solution for different window lengths that 
overlap and differ by one or two frames. 

2. Compare the solution quality, i.e., the value of the 
20 objective function, for two adjacent window lengths. 

3. If there is a predefined improvement in either direction, 
we then, for stability, repeat the exercise for a shorter or 
longer on the first try. If there is consistent information, we 
adjust the window size in the indicated direction. This can be 

25 given a well defined mathematical formula in terms of the 
assignment problems of different dimensions. 
Algorithm Enhancements due to Data Structures (ABP) 
Construction and enumeration of the Best K-Solutions. 
New Nonsmooth Optimization Techniques. 

30 

Sensitivity Analysis 

For sequential processing in which the two dimensional 
assignment algorithm is used, one can use the LP sensitivity 
analysis theory and easily obtain the corresponding answers. 

35 For the multiframe processing, the optimal solution is not 
available; however, there are still two approaches one can 
use. First, the basic algorithm can perform the same sensi- 
tivity analysis at each stage (loop) of the algorithm as is done 
in the two dimensional assignment problem since the evalu- 

4Q ation of function ^.j^ is equivalent to a two dimensional 
assignment problem. Alternately, one could use an LP relax- 
ation of the assignment problem and base the sensitivity on 
the resulting LP problem. We currently see this as an 
important step in finding even better solutions to the assign- 

45 ment problem if so desired. 

New Auction Algorithms (ABP & PS). In this section we 
present a new auction algorithm for the two dimensional 
assignment problem. 

An important step in solving N-dimensional assignment 

50 problems for N^3, is finding the optimal solution of the 
2 -dimensional assignment problem. In particular, we wish to 
solve the 2-dimensional problem which includes the zero- 
index. T3.pically this problem can be thought of as trying to 
find either a minimum cost or maximum utility in assigning 

55 person to objects, tenants to apartments, or teachers to 
classrooms. We will follow the work of Bertsekas and call 
the first index set persons and the second index set objects. 
The statement of this problem is given below (1) when the 
number of persons is m and the number of objects is n. 

60 

m m (1) 

Minimize ££ c ';*'> 

i»0 J'O 
N 

65 Subject to £ xtj B 1 for all J * 1, .... m 
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-continued 



Xij e 10, 1} 



There are a couple of assumptions which we will make 
about (1). First of all* we assume that c l0 and c 0J are 
well-defined and the corresponding variables Xi 0 Xo, fr ee to 10 
be assigned. Second if a cost c jy happens to be undefined- 
then the corresponding variable x,> X to 0. In effect if c,-, is 
undefined, we simply remove this cost and variable from 
inclusion in the problem. 

Notice that there are no constraints on the number of times 15 
person 0 and object 0 can be assigned. But notice that the 
first constraint requires that each non-zero person i must be 
assigned exactly once. Similarly, the second constraint 
forces each non-zero object j to be assigned exactly one 
time. Finally, the last constraint gives an integer solution. 2 o 
although we will see shortly that this constraint can be 
relaxed to admit any solution X^O, One reason for requir- 
ing that all of the costs of the form c ro and c 0J be defined is 
so that we are guaranteed a feasible solution exists for the 
given problem. 25 

Relaxation of One Constraint 

Relaxing the last constraint, we obtain the following 
problem: 



in n — , 1 m 
mm n 



#*0 j*Q 



>0 



Subject to ^ xii m 1 for all i = 1, . . . , m 
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This lower bound is achieved and the second constraint in 
the original problem is satisfied if the following conditions 
are met: 



For i°0, xoj 1 



For i*G, xu B 



1 if coj + u) £0 



0 otherwise 

1 if j e arg min{ett + u\ \ k e 0, 1 n\ 
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(0 otherwise 
All of fs are assigned only one time. 

The relaxation of the first constraint is analogous and would 55 
lead to similar results with i and j exchanged and the 
introduction of the multiplier u\ 
Relaxation of Both Constraints 

This time we will relax both constraints and using duality 
theory obtain a relationship between the multipliers u 1 and ^ 
u 2 . 

a bunch of equations, theorems and the such here. 

From the above relaxation, we arrive at the following 
complementary slackness conditions which have been 
relaxed by a small value e. 65 
Definition: An assignment S and a multiplier set (u\u 2 ) are 
said to satisfy e-complementary slackness (e-CS) if 
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Ctf+uf+u?^ for all (ij)eS, 
c^+i/, for all (ij)eA. 

Forward Auction Algorithm 

(1) Select any unassigned person i 

(2) Determine the following quantities: 

>r=arg min \c A +u k 2 \k«A(t)} 

In the selection of j, above, if a tie occurs between 0 and 
any non-zero index k. then select j, as k. Otherwise, if there 
is a tie between two or more non-zero indices, the choice of 
j, is arbitrary. Also if A(i) consists of only one element, then 
set w=°o. 

(3) Update the multipliers and the assignment: 
If jr=0. then 

(a) Add (i.O) to the assignment. 

(b) Update u/: +-c t0 . 
If jp£(h then 

(a) Add (i.j,) to the assignment 

(b) Remove (i\j f ) from the assignment if j, was previously 
assigned. 

(c) Update u/: ^-Kw.-v^sw.-c^+e. 

(d) Update u/: ^(c^+u^^w^e. 
Reverse Auction Algorithm 

(1) Select any unassigned object j. such that Cq^+u/^. 

(2) Determine the following quantities: 

ij=arg mm {c Ji +u k l \k€B0')} 

In the selection of i^. above, if a tie occurs between 0 and 
any non-zero index k. then select j ( as k. Otherwise, if there 
is a tie between two or more non-zero indices, the choice of 
j, is arbitrary. Also if B(j) consists of only one element, then 
set jj 00. 

(3) Update the multipliers and the assignment: 
If i=0, then 

(a) Add (O.j) to the assignment. 

(b) Update u 2 : =>-c 0/ 
If i/*0. then 

(a) Add (i y .j) to the assignment 

(b) Remove (i^j') from the assignment if (i^.j) was previ- 
ously assigned. 

(c) Update u,/: =u l /-+<Y r P y >K=^-c iy +e 

(d) Update u/: ^cyHi.^y^ 
Combined Forward/Reverse Auction Algorithm 

1. Assume that u 2 is given as an arbitrary multiplier. 

2. Adjust the value of u 2 for each object j as follows: 



If Coyfu/O, then set u, 2 : 



3. Run iterations of the Forward Auction Algorithm until 
all persons become assigned. 

4. Run iterations of the Reverse Auction Algorithm until 
all of the objects become assigned. 

Note at the completion of the Forward auction step we 
have the following conditions satisfied: 



